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Summary 

This  study  had  two  goals,  the  first  of  which  was  to  examine  approaches  for  calculating  the  hydrodynamics 
of  breaking  waves  and  to  assess  the  impact  of  hydrodynamic  model  errors  on  the  prediction  of  radar 
backscatter.  To  address  this,  Reynolds-averaged  Navier-Stokes  (RANS)  computations  of  stationary 
hydrofoil-generated  breaking  waves  were  carried  out,  including  the  modeling  of  the  breaking  region. 
These  results  were  compared  to  experimental  data.  To  determine  the  impact  of  the  errors  in 
hydrodynamic  modeling  on  predictions  of  radar  backscatter,  a  subset  of  the  hydrodynamic  results  were 
used  as  input  to  the  Veridian  Scattering  Model  (VSM)  and  the  results  were  compared  to  those  of  Ericson 
et  al.  (1999)  which  uses  similar  backscatter  modeling  but  with  experimentally  measured  hydrodynamic 
inputs.  The  second  goal  of  this  study  was  to  use  the  information  gained  through  this  effort  to  define  the 
research  needs  in  this  area. 

The  conclusions  of  the  study  are  that  computations  of  breaking  waves  produce  qualitatively  correct 
behavior,  in  the  sense  that  energy  is  dissipated,  turbulence  is  produced  in  the  breaking  region  and  the 
wave  amplitude  is  reduced;  however,  the  results  are  not  particularly  accurate  and  the  present  methods  are 
not  robust.  For  some  cases  the  wave  amplitude  was  over-estimated  and  for  some  it  was  under-estimated, 
indicating  incorrect  amounts  of  dissipation  from  the  breaking  model,  while  for  some  cases,  a  converged 
solution  could  not  be  obtained  at  all.  For  a  subset  of  the  cases,  the  VSM  radar  backscatter  model  was 
used  to  estimate  both  the  RMS  surface  fluctuation  level  and  the  radar  cross-section  (RCS)  for  X-Band  at 
45  degrees  incidence.  The  RMS  surface  fluctuations  in  all  cases  drop  too  slowly  as  the  short  waves 
propagate  downstream  from  the  breaking  region.  This  is  traced  to  the  fact  that  the  RANS-predicted 
velocity  field  in  the  breaking  region  is  not  accurate.  The  radar  cross-section  was  estimated  using  the 
small-perturbation  method  (SPM)  and  the  results  were  compared  to  those  generated  using  the  SPM  model 
with  the  correct  hydrodynamic  inputs,  derived  from  experimental  observations,  presented  in  Ericson  et  al. 
(1999).  Overall,  the  peak  RCS  values  for  the  breaking  region  tend  to  be  over-estimated  by  as  much  as 
12  dB,  and  the  downstream  RCS  values  tend  to  be  over-estimated  by  as  much  as  25  dB  when  compared  to 
the  results  of  the  SPM  model  using  the  correct  hydrodynamic  inputs.  Hence,  hydrodynamic  errors 
incurred  in  modeling  breaking  waves  can  lead  to  significant  errors  in  radar  backscatter  estimates. 

A  number  of  areas  where  research  is  needed  have  been  identified.  In  order  to  predict  the  behavior  of 
breaking  waves  for  flows  of  practical  interest,  approaches  such  as  the  Reynolds-averaged  Navier-Stokes 
equations  must  be  employed.  This  type  of  approach  requires  that  some  of  the  details  of  the  flow  must  be 
represented  by  models  which  need  to  accurately  reflect  the  important  physical  processes.  The  results  of 
this  study  indicate  that  improvements  must  be  made  in  virtually  all  aspects  of  modeling  wave  breaking. 
These  include  modeling  the  onset  of  wave  breaking,  modeling  the  generation  of  surface  disturbances  in 
the  breaking  region,  predicting  the  propagation  of  short  waves  through  and  beyond  the  breaking  region, 
and  development  of  proper  free-surface  boundary  conditions.  Another  important  research  area  is 
development  of  more  robust  and  efficient  computational  approaches,  including  radiation  boundary 
conditions  and  efficient  solvers. 

To  undertake  the  development  of  improved  modeling  approaches,  new,  more  detailed  information  on  the 
behavior  of  breaking  waves  is  required.  There  are  two  complementary  ways  in  which  this  information 
can  be  obtained:  The  first  is  through  direct  numerical  simulations  (DNS)  based  on  numerical  solution  of 
the  exact  governing  equations  for  the  three-dimensional  time-evolution  of  breaking  waves.  DNS  allows 
the  details  of  the  waves  to  be  examined  both  qualitatively  and  quantitatively,  and  can  provide  access  to 
un-measurable  quantities  such  as  the  subsurface  pressure  fluctuations.  These  detailed  data  can  be  used  to 
guide  the  development  of  the  needed  models.  The  drawback  to  this  approach  is  that  presently  it  can  only 
be  applied  to  relatively  short  breaking  waves  (i.e.  waves  substantially  less  than  a  meter  in  wavelength) 
and  even  then  it  taxes  available  computing  resources.  The  second  approach  for  obtaining  the  needed 
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information  is  experimental  studies  of  breaking  waves.  Here,  a  more  limited  set  of  information  can  be 
obtained  using  available  measurement  technology,  but  this  can  be  done  for  the  range  of  scales  from  DNS 
to  scales  approaching  situations  of  practical  interest.  Detailed  measurement  of  surface  structure,  velocity 
and  subsurface  turbulence  are  needed.  Through  a  coordinated,  combined  approach  of  DNS  and 
experiments,  the  information  necessary  to  develop  first-principles  models  for  wave  breaking  and  the 
associated  small-scale  disturbances  important  for  radar  backscatter  can  be  obtained. 

A  final  issue  is  the  framework  within  which  modeling  of  breaking  waves  should  be  addressed.  Reynolds- 
averaged  approaches  do  not  exclude  any  of  the  physical  processes  which  naturally  occur  in  wave 
breaking,  hence  RANS  approaches  should  be  able  to  predict  wave  breaking  accurately,  given  appropriate 
models.  For  robustness,  it  appears  that  a  computational  approach  which  admits  overturning  waves 
(multiple-valued  free  surface),  such  as  a  level-set  or  an  improved  surface-tracking  approach,  will  be 
needed.  Potential  flow  methods  impose  some  requirements  on  the  flow  (inviscid,  irrotational)  which  are 
violated  in  breaking  waves;  however,  potential  flow  methods  are  significantly  less  demanding 
computationally.  For  this  reason,  potential  flow  methods  would  be  attractive,  if  they  could  be  extended  to 
model  wave  breaking  without  resorting  to  ad-hoc  approaches.  This  should  be  investigated  as  well. 


2 


DW200392.8.F.Nov04 


Introduction 


Breaking  waves  on  the  ocean  surface  can  yield  a  significant  radar  return;  however,  the  prediction  of  radar 
backscatter  from  breaking  waves  is  limited  by  our  ability  to  predict  the  details  of  the  mean  velocity  field 
and  surface  shape  as  well  as  the  small-scale  hydrodynamic  disturbances  which  result  from  wave  breaking. 
The  first  goal  of  this  study  was  to  examine  approaches  for  calculating  the  hydrodynamics  of  breaking 
waves  and  to  assess  the  impact  of  hydrodynamic  model  errors  on  the  prediction  of  radar  backscatter.  To 
address  this,  first  the  hydrodynamic  models  were  used  to  calculate  the  mean  surface  elevation  and 
velocity  fields  for  two-dimensional  stationary  hydrofoil-generated  waves.  These  results  were  compared 
to  available  experimental  data.  To  determine  the  impact  of  the  errors  in  hydrodynamic  modeling  on 
predictions  of  radar  backscatter,  a  subset  of  the  hydrodynamic  results  were  used  as  input  to  the  Veridian 
Scattering  Model  (VSM).  The  VSM  radar-cross-section  results  could  then  be  compared  to  data  to  assess 
the  end-to-end  performance  of  the  combined  hydrodynamic-radar-backscatter  modeling  approach. 

The  second  goal  of  this  study  was  to  use  the  information  gained  through  this  effort  to  define  the  research 
needs  in  this  area.  To  calculate  radar  backscatter,  the  two-dimensional  energy  spectrum  of  the  small-scale 
waves  is  needed,  along  with  the  mean  surface  slope.  The  ability  to  calculate  these  quantities  will  be 
assessed,  and  the  deficiencies  in  the  approaches  will  be  identified.  This  will  identify  areas  where  new 
approaches  or  new  knowledge  are  required. 

Hydrodynamic  predictions  of  waves  can  be  made  via  numerical  solution  of  the  Navier-Stokes  equations 
subject  to  the  appropriate  free-surface  boundary  conditions.  For  practical  flows  which  are  steady  in  the 
mean,  simplifications  to  these  equations  must  be  made  in  order  to  make  the  computations  feasible.  Two 
approaches  for  this  exist:  potential  flow,  where  viscous  and  rotational  effects  are  formally  eliminated 
from  the  equations,  and  Reynolds-averaged  approaches,  where  the  equations  are  time-averaged.  Potential 
flow  calculations  are  significantly  faster,  but  cannot  account  for  viscous  effects  or  turbulence.  Reynolds- 
averaged  Navier-Stokes  (RANS)  calculations  require  significantly  more  computational  effort,  as  well  as 
turbulence  modeling,  but  are  a  framework  which  can  accommodate  both  viscous  and  turbulence  effects. 

In  typical  calculations  of  waves  with  either  RANS  or  potential  flow,  wave  breaking  is  not  considered  and 
solutions  can  be  obtained,  even  for  very  steep  waves  which  in  the  real  world  would  break.  In  this  case, 
the  surface  elevation  and  velocity  fields  predicted  are  solutions  to  the  equations,  but  are  incorrect,  and 
radar  backscatter  predictions  based  on  these  results  will  also  be  incorrect.  If  the  waves  become  steep 
enough,  either  RANS  or  potential  flow  computations  will  fail  unless  a  breaking  model  is  invoked. 
Modeling  the  effects  of  wave  breaking  is  presently  relatively  crude,  relying  on  corrections  to  the  surface 
pressure  and  velocity  that  stabilize  the  steep  wave  and  prevent  overturning. 

The  small-scale  unsteady  surface  disturbances  generated  by  wave  breaking  are  important  for  radar 
backscatter;  however,  modeling  the  generation  of  short  waves  by  wave  breaking  has  not  been  addressed 
in  a  significant  way  by  the  hydrodynamics  community.  An  empirical  model  for  the  short-wave  spectrum 
in  the  breaking  region  of  a  stationary  wave  has  been  developed  (see  e.g.  Walker  et  al.  1996).  This 
empirical  spectrum  is  used  in  VSM  (Ericson  &  Wackerman  2000)  to  prescribe  the  behavior  of  the  short 
waves  in  the  breaking  region.  The  propagation  of  short  waves  as  they  leave  the  breaking  region  and 
interact  with  the  surrounding  surface  velocity  field  is  calculated  via  solution  of  the  wave-action  balance 
equations  (Komen  et  al.  1994),  which  then  provides  the  short-wave  spectrum  for  the  entire  surface.  To 
accomplish  this,  VSM  requires  a  description  of  the  large-scale  mean  hydrodynamics  (surface  elevation, 
surface  velocity)  along  with  a  map  of  the  breaking  regions.  These  inputs  are  derived  from  either  RANS 
or  potential-flow  solutions.  An  alternative  approach  is  to  specify  a  turbulence-related  short-wave  source 
term  directly  in  the  action  balance  equations  (Walker  2000).  This  has  not  been  applied  to  breaking 
waves,  but  has  been  used  in  other  flows. 
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Given  the  short-wave  spectrum,  the  radar  backscatter  can  be  estimated  using  different  approximate 
theories.  The  comparisons  to  be  made  in  this  study  will  be  for  laboratory  measurements  at  moderate 
incidence  angles  (±45°).  For  this  case  the  small-perturbation  method  (SPM),  as  described  by  Ericson  et  al 
(1999)  is  used. 

In  the  following  sections,  first  the  hydrodynamic  modeling  is  described  and  issues  of  model  type  and  grid 
resolution  are  addressed.  Then  hydrodynamic  results  calculated  for  breaking  waves  are  shown  and 
compared  to  available  experimental  data.  The  hydrodynamic  results  are  then  used  as  a  basis  for  radar 
backscatter  predictions  and  again  compared  to  experimental  data. 
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Approach 


Hydrodynamic  model  calculations  were  carried  out  mainly  using  a  Reynolds-averaged  Navier-Stokes 
(RANS)  approach;  however,  as  a  check  on  the  accuracy  of  those  solutions,  potential  flow  calculations 
were  performed  as  well  for  some  cases.  Results  from  two  different  RANS  codes,  CFDSHIP-IOWA  and 
UNCLE,  which  have  different  numerical  formulations,  were  used,  so  as  to  rule  out  systematic  errors.  The 
RANS  solutions  for  some  of  the  cases  examined  diverged,  because  the  predicted  waves  became  too  steep 
and  overturned.  To  deal  with  these  cases,  a  wave-breaking  model,  which  redistributes  momentum  and 
dissipates  energy  in  the  waves  was  used.  Finally  the  RANS  results  for  a  subset  of  the  cases  were  used  as 
input  to  radar  backscatter  computations  using  VSM.  The  approaches  used  for  each  of  these  are  described 
here. 


RANS  Computations 

Reynolds  Averaged  Navier-Stokes  (RANS)  governing  equations  are: 

Vu  =  0  (1) 


p| - l-u-Vu  =  -Vp  +  pV2u-  V-t,  +pg 

Ot  ) 


(2) 


where  x,  is  the  Reynolds  stress.  The  Reynolds  stress  is  calculated  either  using  the  k-e  model  from  the 
Boussinesq  hypothesis  or  directly  from  a  Reynolds-stress  model  (RSM).  Boundary  conditions  imposed  on 
the  free  surface  is  the  kinematic  boundary  condition 

— “  +  u*Vr|  =  w  on  z  =r|(x,  y,t) ,  (3) 

dt 

and  the  dynamic  boundary  condition 


-  pn  +  \i(Vu  +  VuT)-n  =  -pAmn.  (4) 

The  solution  of  the  RANS  equations  is  obtained  using  two  Computational  Fluid  Dynamics  (CFD)  codes: 
CFDSHIP-IOWA  and  UNCLE. 

CFDSHIP-IOWA  is  based  on  the  Finite  Difference  Method  (FDM)  that  solves  the  equations  in  three- 
dimensional  space  and  unsteady  mode  (Paterson  et  al.,  2003).  It  includes  a  free  surface  tracking 
capability,  where  the  computational  grid  is  fitted  to  the  free-surface  boundary.  Turbulence  is  modeled 
with  a  two-equation  k-e  model.  The  code  is  being  used  primarily  for  the  study  of  ship  hydrodynamics  and 
is  considered  a  state-of-the-art  code  for  this  type  of  problems.  The  code  is  parallelized  and  can  be  run 
using  Message  Passing  Interface  (MPI). 

UNCLE  is  based  on  the  Finite  Volume  Method  (FVM).  The  equations  are  solved  in  three-dimensions 
and  the  simulations  can  be  performed  both  in  steady  and  unsteady  modes.  The  convection  terms  are 
discretized  using  the  Roe’s  flux-splitting  approach,  and  the  momentum  equations  are  solved  in  a  coupled 
manner.  Turbulence  is  modeled  either  with  the  k-e  or  full  Reynolds  Stress  Model  (RSM).  The 
calculation  of  the  turbulence  model  is  decoupled  from  the  calculation  of  the  velocity  components.  It  also 
has  the  surface  tracking  capability  and  can  be  run  in  parallel  mode  using  MPI. 
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Potential  Flow  Computations 


In  case  of  incompressible,  inviscid,  irrotational  flow,  it  is  possible  to  define  a  potential,  O,  so  that  the 
velocity  is  given  by  the  gradient  of  a  scalar  potential  function, 

u  =  VO. 

Using  such  a  definition  of  the  potential,  the  governing  equations  can  be  rewritten  as: 

V2O  =  0 

30 

—  +  p  +  ^pVO- VO  +  pgz  =  C(0 

3 1 

Consequently,  the  kinematic  and  dynamics  boundary  conditions  become  equal  to: 

— -  +  V0-Vt|  =  0  on  z  =  t|(jc,  y,t) 


The  potential  flow  equations  are  solved  using  a  de-singularized  potential  method  (Scullen  &  Tuck,  1995; 
Scullen,  1998).  The  effects  of  the  submerged  hydrofoil  are  implemented  with  a  linearly  varying  vortex 
panel  method  (Kuethe  &  Chow,  1986).  This  implementation  is  described  in  detail  in  the  Appendix, 
below. 


Trough  B 


Looking 

Upwave 


i  Lc 

Downwavl 


•1 


BWiS 


Breaking  Model 


For  hydrofoil-generated  waves,  a  trough  forms  above  the  hydrofoil,  followed  by  a  wave  train.  If  the 
disturbance  caused  by  the  foil  is  large  enough,  due  either  to  the  depth  of  the  foil  or  its  angle  of  attack,  the 
initial  wave  crest  will  break.  Typical  pressure  and  velocity  fields  in  such  a  case  are  shown  in  Figure  la) 
and  Figure  lb).  Theory  indicates  that  the  orbital  velocities  associated  with  the  waves  increase  in 
proportion  to  the  wave  amplitude.  When  the  orbital  velocity  at  the  wave  crest  exceeds  the  phase  velocity 
of  the  waves,  the  waves  will  break.  The  process  can  be  very  abrupt  as  in  the  case  of  a  plunging  breaker, 
or  more  tempered  as  a  spilling  breaker  (Figure  lc  and  Figure  Id).  The  spilling  breaker  causes  both 
velocity  and  pressure  modification  close  to  the  wave  surface.  Simulation  of  the  wave  under  breaking 
conditions  requires  the  inclusion  of  these  effects. 


Figure  2.  Schematic  of  the  breaker  and  the  definitions  of  the  modeling  parameters. 


A  steady  breaker  can  be  modeled  using  the  semi-analytical  approach  of  Cointe  &  Tulin  (1994).  In  this 
approach,  the  breaking  region  (the  region  of  partially  aerated  white  water  forward  of  the  breaking  wave 
crest)  is  assumed  to  ‘ride’  on  of  the  forward  face  of  the  wave.  At  the  forward  edge  of  the  breaking  region, 
the  ‘toe’  of  the  breaker,  the  free-surface  surface  streamline  becomes  submerged  under  the  breaking  region 
and  re-emerges  at  the  wave  crest  (see  the  blue  line  in  Figure  2).  The  weight  of  the  breaking  region 
increases  the  pressure  along  the  free-surface  streamline.  If  one  assumes  a  shape  and  density  for  the 
breaking  region,  and  further  that  the  total  head  is  conserved  along  the  streamline,  the  velocity  and 
pressure  along  the  free-surface  streamline  under  the  breaking  region  can  be  determined.  These  quantities 
are  then  used  in  place  of  the  usual  free  surface  kinematic  and  dynamic  boundary  conditions,  (3)  and  (4), 
respectively.  From  a  momentum  conservation  standpoint,  the  force  associated  with  the  increased  pressure 
acting  on  the  free-surface  streamline  results  in  an  increased  drag  on  the  hydrofoil,  due  solely  to  wave 
breaking.  The  approach  as  implemented  can  be  described  in  reference  to  the  schematic  of  the  breaking 
wave  crest,  shown  in  Figure  2.  Following  the  Rhee  &  Stem  (2000)  implementation,  the  breaking  model 

is  invoked  for  any  wave  crests  which  exhibit  0(x)  >15°  anywhere  between  the  trough  and  crest.  Here, 
the  blue  line  represents  the  free-surface  streamline  5.  The  dotted  line  in  Figure  2  represents  the  observed 
geometry  of  the  breaking  region  which  is  approximated  by  a  triangle  (the  red  line).  The  toe  of  the  breaker 
is  located  at  the  point  labeled  a,  and  its  location  can  be  specified  in  different  ways;  here,  we  will  follow 
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Muscari  &  Di  Mascio  (2003)  and  define  the  vertical  distance  between  the  crest  and  the  toe  (points  b  and 
a)  as 


h*=zb-za=0MH, 


(10) 


where  the  choice  of  the  coefficient  0.64  by  Muscari  &  Di  Mascio  was  based  on  the  experimental  data  of 
Duncan  (1983).  It  is  assumed  that  there  is  a  discontinuity  in  the  velocity  at  a,  given  by 

=  (11) 

where  (3  is  a  ‘mixing’  parameter  with  values  between  0.5  and  0.7,  and  the  value  of  u0  ( xa )  is  determined 
by  applying  Bernoulli’s  equation  on  the  streamline  s  upstream  of  a.  This  yields 

«0  (xa)=c2 -2gz{xa),  (12) 


where  c  is  the  phase  velocity  of  the  wave  (i.e.  the  flow  velocity  far  upstream  of  the  foil).  In  the  breaking 
region,  the  pressure  on  the  surface  between  a  and  b  is  assumed  hydrostatic 


P(x )  =  p g  I  z(x)  +  h(x) 


(13) 


The  ratio  pBlv  /p  represents  the  ratio  between  the  density  of  the  water/air  mixture  within  the  breaker  and 
pure  water  density.  The  value  for  this  ratio  was  estimated  by  Cointe  &  Tulin  (1994)  to  be 
0.2  —  p  gw  /  p  —  0.6  where  the  lower  limit  represents  incipient  breaking  and  the  top  limit  the  strong 
breaking  process;  a  value  of  0.5  was  used  in  the  calculations.  The  height  of  the  triangular  breaker  is 
calculated  as 


h(x)  =  z{xb)-z(x).  (14) 

Streamline  velocity  for  any  point  between  a  and  b  is  calculated  using  the  Bernoulli  equation  and 
equations  (1 1)— (14): 


u]  (*)  =  u]  (xa )  +  2  [p(xa )  -  p(xj\ .  (15) 

By  using  the  local  surface  slope,  6,  and  the  streamline  velocity  us ,  the  two  surface  velocity  components 
are  calculated  as: 


u(x)  =  us  (x)  cos  0(jc)  ,  w(jc)  =  u s  (x)  sin  9(x)  (16) 

The  surface  pressure  from  (13)  and  the  surface- velocity  components  (15)  are  then  used  to  specify  the 
boundary  conditions  in  the  breaking  region. 

Short-Wave  Modeling 

Reynolds-averaged  or  potential-flow  approaches  which  can  be  used  to  calculate  the  mean  surface- 
elevation  and  velocity  fields;  however,  the  short  waves  which  are  responsible  for  radar  backscatter  are 
not  resolved  using  these  approaches.  The  short-wave  spectrum  is  modeled  using  VSM  (Ericson  & 
Wackerman  1999).  In  that  approach,  the  breaking  regions  are  identified  using  a  vertical  acceleration 
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criterion,  an  empirical  short-wave  spectrum  is  specified  in  the  breaking  region,  and  the  evolution  of  the 
short-wave  spectrum  outside  the  breaking  region  is  calculated  using  the  wave-action  balance  equations. 
The  breaking  criterion  is 


^<0.15g  (17) 

which  identifies  regions  of  large  downward  acceleration.  In  these  regions  the  short-wave  spectrum  is 
assumed  isotropic  and  is  set  to 


5(k)  =  0.004ork-35  k>2x/Lb 

(18) 

=  0  k<2x/Lb 

where  k  is  the  wavenumber  and  Lb  is  the  length  of  the  breaking  region  in  the  wave-propagation  direction. 
The  factor  a  is  given  by 


a = exp 


3.5  (  2#  v 
4 


(19) 


and  is  included  to  yield  zero  slope  in  the  spectrum  at  k  =  2x1  Lh .  It  only  affects  the  spectrum  very  near 
the  cutoff.  The  propagation  of  the  waves  outside  the  breaking  region  is  done  using  the  wave  action- 
balance  equations  which  govern  the  evolution  of  the  short-wave  spectrum 

za  c 

Y+V'(C  A)  =  -  (20) 

dt  g  <3 

where  A  is  the  wave-action  spectrum  (the  energy  spectral  density  divided  by  the  wave  frequency),  Cg  is 

the  vector  group  velocity  (the  propagation  velocity  for  the  wave  energy  in  both  physical  and  spectral 
space)  and  the  source  term  S  represents  the  effects  of  wind  growth,  viscous  dissipation,  etc. 

The  identification  of  the  breaking  region  is  a  key  element  of  the  prediction  of  the  short-wave  spectrum,  as 
it  defines  both  the  location  of  the  breaking  regions  and  the  length  of  the  breaking  region,  Lb .  For  the 
radar-cross-section  comparisons  below,  the  breaking  region  was  determined  in  two  ways.  One  was  using 
the  native  VSM  approach  (based  on  the  acceleration).  The  second  approach  used  the  information  from 
the  breaking  model  implemented  in  the  RANS  code  (based  on  slope). 

Radar  Backscatter  Model 

The  radar  backscatter  model  used  for  this  study  was  the  small-perturbation  method  (SPM),  described  in 
Ericson  et  al.  (1999).  The  SPM  solution  for  radar  backscatter  is  given  by  Valenzuela  (1978)  as 

Op 7  =87tko4|Yp/,(0,)|2['P(k6)  +  'P(-kfc)],  (21) 

where  T  is  the  surface-elevation  energy  spectral  density,  ypp  is  a  polarization-dependent  reflection 

coefficient,  0(.  is  the  local  incidence  angle,  k0  is  the  radar  wavenumber,  and  kb  =  2 k0  sin  0,  is  the  Bragg 
wavenumber. 
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RANS  Implementation 

In  this  section  the  implementation  issues  unique  to  the  hydrofoil-generated  waves  examined  here  are 
discussed.  First  wave  reflections  are  discussed,  along  with  methods  for  eliminating  them.  Then  the 
computational  domain  is  described.  Finally  the  grid  requirements  are  examined. 

The  results  indicate  that  in  order  to  obtain  a  RANS  solution,  the  domain  must  be  large  and  the  grid 
sufficiently  stretched  eliminate  reflection  from  the  upstream  and  downstream  domain  boundaries. 
Additional  damping  of  the  waves  is  sometimes  needed.  What  was  also  learned  is  that  as  long  as  the  flow 
about  the  foil  is  well  resolved,  the  near-field  solution  for  the  free  surface  elevation,  and  hence  the  velocity 
field,  will  be  accurately  captured.  If  the  grid  resolution  away  from  the  foil  is  not  adequate,  the  far  field 
waves  will  be  attenuated.  This  near-field  accuracy  vs.  far-field  attenuation  has  been  seen  in  Kelvin-wake 
predictions,  and  is  tied  directly  to  grid  resolution  and  its  impact  on  wave  propagation. 

Wave  Reflections  from  Domain  Boundaries 


The  major  issue  in  the  simulation  proved  to  be  the  reflection  of  waves  from  both  upstream  and 
downstream  boundaries.  The  flow  domain  is  initialized  with  the  flat  surface.  Waves  start  to  form  during 
the  initial  transient  and  the  wave  train  behind  the  hydrofoil  moves  towards  the  downstream  boundary. 
Similarly,  a  single  wave  is  generated  continuously  in  front  of  the  foil  and  propagates  towards  the 
upstream  boundary.  Once  the  waves  reach  the  boundaries,  they  reflect  and  cause  considerable 
oscillations  in  wave  amplitudes  within  the  domain.  These  oscillations  are  further  enhanced  at  the 
boundaries  and  cause  severe  mass  loss.  Despite  the  fact  that  the  calculations  can  continue  in  this  mode, 
the  convergence  rate  is  reduced  and  the  results  lose  physical  meaning. 

Several  authors  observed  the  same  issue  in  simulating  the  free  surface  flows  (Iafrati  et  al.,  2001; 
Brummelen  et  al.,  2001;  Muscari  &  Di  Mascio,  2003).  They  proposed  mesh  extensions  away  from  the 
region  of  interest  thus  resulting  in  increased  numerical  dissipation.  In  addition,  both  Iafrati  et  al.  (2001) 
and  Muscari  &  Di  Mascio  (2003)  proposed  wave  damping.  We  used  both  methods  to  prevent  the 
reflected  waves.  Mesh  size  was  increased  geometrically  away  from  the  region  of  interest  resulting  in 
cells  a  few  hydrofoil  lengths  in  size.  Wave  damping  was  implemented  as  well: 


dn  dn  d7j 
—L+u-L  +  v-*- 
ot  dx  dy 


w-drj 


(22) 


where  is  7]  surface  elevation  and  d  is  the  damping  factor.  The  damping  factor  was  set  to  zero  in  the  region 
of  interest  and  increased  geometrically  away  from  it. 

Computational  Domain 

Figure  3  shows  details  the  computational  grid  or  mesh.  The  mesh  had  increasingly  larger  cells  away  from 
the  hydrofoil  towards  both  ends  of  the  computational  domain,  as  shown  in  Figure  3a)  and  Figure  3c).  The 
most  commonly  used  mesh  consisted  of  four  blocks,  each  with  400x20  nodes  in  the  plane.  As  can  be 
seen  in  Figure  3b)  the  mesh  around  the  hydrofoil  was  non-uniform  with  refinement  in  the  boundary  layer 
and  at  the  front  and  trailing  edges  of  the  hydrofoil.  Flow  was  quasi  two-dimensional  with  symmetry 
boundary  conditions  imposed  in  the  ^-direction.  Even  though  the  considered  geometry  was  two- 
dimensional,  both  UNCLE  and  CFDSHIP-IOWA  can  perform  only  three-dimensional  calculations.  As  a 
result,  simulations  performed  with  UNCLE  require  meshes  at  least  3  cells  thick,  and  CFDSHIP-IOWA  at 
least  5  cells  thick. 
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Figure  3.  Mesh  used  in  the  simulation:  a)  mesh  extension  before  the  hydrofoil,  b)  mesh  detail  around  the 
hydrofoil,  c)  mesh  extension  after  the  hydrofoil. 


The  entire  wave  profile  calculated  for  the  non-breaking  wave  case  measured  by  Duncan  (1983)  is  shown 
in  Figure  4.  The  hydrofoil  submergence  depth  in  this  case  was  0.261  m.  Observed  can  be  that  the  surface 
is  flat  on  both  ends  of  the  computational  domain.  This  is  a  result  of  using  wave  damping  and  large  cell 
sizes  that  induce  large  numerical  diffusion. 
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Figure  4.  The  calculated  wave  profile  for  the  Duncan  0.261  m  case  shown  for  the  entire  computational 
domain. 


- T— 

. . 1 

II 

ls}«  --- 

I 

ZZJ 

—  Duncan  (1983);  /i=0.261  m 
--  Calculation  j 

Grid  Variations 

In  order  to  establish  the  accuracy  of  the  computations  to  be  presented  below  for  breaking  waves,  the  CFD 
codes  were  applied  to  non-breaking  wave  cases.  In  these  flows,  the  grid  requirements  and  overall 
behavior  of  the  computations  were  assessed,  independent  of  the  complications  associated  with  modeling 
wave  breaking.  Simulations  of  the  non-breaking  waves  were  performed  using  the  two  RANS  codes, 
CFDSHEP-IOWA  and  UNCLE.  First,  the  sensitivity  of  the  RANS  simulations  to  the  computational  mesh 
was  investigated.  Next,  the  description  of  the  RANS  calculations  and  the  details  of  the  computational 
domain  and  other  issues  addressed  in  carrying  out  the  RANS  computations  are  presented. 
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It  should  be  noted  that  the  grid  studies  described  in  this  section  were  carried  out  using  a  slightly  imperfect 
geometry  for  the  NACA  0012  hydrofoil  used  in  the  experimental  work.  This  has  no  bearing  on  the  results 
of  the  grid  study.  All  the  other  results  were  for  the  precise  geometries  used  in  the  experiments. 

Two  mesh  types  were  used  in  the  RANS  simulations  to  verify  calculation  sensitivity  to  the  mesh  type.  A 
number  of  Cartesian  meshes  were  used  to  produce  most  of  the  results  generated  for  this  study.  An 
example  of  a  Cartesian  mesh,  shown  in  Figure  5a)  had  fine  grid  spacing  in  the  boundary  layer,  and  at  the 
leading  and  trailing  edges  of  the  hydrofoil.  This  Cartesian  mesh  had  32,000  points  in  the  x-z  plane.  In 
addition  to  using  different  grid  densities  for  the  Cartesian  meshes,  a  single  C-mesh,  shown  in  Figure  5b), 
was  used.  In  the  C-mesh,  the  grid  wraps  around  the  foil  and  is  generally  more  orthogonal.  The  C-mesh 
was  refined  near  the  hydrofoil  boundary,  and  it  had  85,364  cells  in  two-dimensions. 


Figure  5.  Meshes  used  in  the  simulations:  a)  Cartesian  mesh  and  b)  C-mesh. 


The  comparison  between  the  calculated  wave  profiles  obtained  on  the  two  mesh  types  together  with  the 
experimental  data  are  shown  in  Figure  6a).  The  predicted  location  and  depth  of  the  trough  are  in  a  very 
good  agreement.  The  result  obtained  on  the  C-mesh  has  a  slightly  higher  crest.  This  is  primarily  due  to  the 
higher  mesh  resolution  in  this  case.  The  predicted  pressure  coefficient,  Cp ,  distributions  along  the 

hydrofoil  are  shown  in  Figure  6b)  Unfortunately,  no  experimental  data  are  available  for  the  pressure 
distribution.  Flow  acceleration  is  higher  along  the  top  hydrofoil  surface  in  the  case  with  the  C-mesh.  As  a 
result,  pressure  drop  along  the  same  surface  is  more  than  for  the  Cartesian  mesh.  Nevertheless,  the 
differences  are  not  large  and  are  not  affecting  the  calculated  wave  profiles  in  a  significant  way.  Due  to  the 
large  number  of  cases  that  had  to  be  calculated,  increase  in  the  computational  time  required  for  the  C- 
mesh  could  not  be  afforded  and  all  the  simulations  were  performed  on  the  Cartesian  meshes. 

The  results  of  a  resolution  study  for  the  Cartesian  meshes  are  shown  in  Figure  7.  Several  calculations 
were  performed  with  both  UNCLE  and  CFDSHIP-IOWA  codes  on  meshes  ranging  from  32,000  grid 
points  in  the  plane  to  121,000  grid  points  in  the  plane.  The  comparison  for  the  wave  profiles  is  shown  in 
Figure  7a).  The  results  are  all  consistent  for  the  initial  trough  depth  and  the  height  of  the  first  crest.  For 
the  coarser  meshes,  the  amplitude  of  the  waves  attenuates  with  downstream  distance.  The  calculated 
results  under-predict  the  depth  of  the  trough  and  the  height  of  the  measured  crest  for  all  cases,  when 
compared  to  the  experimental  data.  The  computational  results  are  almost  identical  inside  the  region 
where  the  measurement  data  is  available.  Similarly,  the  predicted  pressure  distributions  on  the  hydrofoil 
surface  shown  in  Figure  7b)  differ  only  slightly.  The  behavior  of  these  results  are  consistent  with  those 
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discussed  by  Reed  &  Milgram  (2000)  for  RANS  calculations  of  Kelvin  wave  patterns.  They  showed  that 
even  though  RANS  computations  could  accurately  reproduce  the  near-field  wave  pattern,  the  waves 
attenuated  rapidly  far  from  the  ship.  The  results  shown  above  indicate  that  this  can  be  traced  to 
inadequate  grid  resolution  resulting  in  far-field  attenuation. 


x  [m]  x  [m] 


Figure  6.  Results  calculated  on  the  Cartesian  and  C-meshes  for  the  Duncan  0.210  m  case  obtained  with 
UNCLE.  Shown  are  a)  wave  profiles  and  b)  surface  pressure  distributions. 
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Figure  7.  Results  calculated  with  CFDSHIP-IOWA  and  UNCLE  for  the  Duncan  0.261  m  case.  Shown 
are  a)  wave  profiles  and  b)  surface  pressure  distributions. 
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Results 


In  this  section  the  results  from  RANS  and  potential-flow  computations  for  non-breaking  are  presented. 
This  is  followed  by  results  RANS  computations  for  breaking  waves  using  the  Cointe  &  Tulin  (1994) 
breaker  model.  The  section  closes  with  radar  cross  section  computations  using  the  RANS  results  as  input. 

For  the  non-breaking  waves,  the  main  results  indicate  that  while  RANS  and  potential  flow  approaches  are 
capable  of  accurately  calculating  the  wavelength  of  hydrofoil-generated  stationary  breakers,  potential 
flow  methods  tend  to  over-estimate  the  wave  amplitude,  while  RANS  appears  to  under-estimate  the 
amplitude.  It  should  be  noted  that  the  amount  of  data  available  for  comparison  is  limited  in  both  extent, 
and  type.  Velocity  fields  along  with  surface  elevation  profiles  would  allow  more  detailed  comparisons 
and  more  complete  validation  of  calculations. 

For  the  breaking  wave  cases,  three  different  data  sets  were  examined  and  the  model  was  shown  to  have 
qualitatively  correct  behavior;  i.e.  the  energy  is  dissipated  and  turbulence  is  produced  in  the  breaking 
region  and  the  wave  amplitude  is  reduced.  The  cases  of  Duncan  (1983)  were  not  successfully  calculated, 
either  with  or  without  the  breaking  model.  These  breakers,  including  the  incipient  breaking  case  could 
not  be  converged  to  a  steady  solution.  Given  that  steady  solutions  are  obtained  for  non-breaking  wave 
cases,  this  points  to  the  breaking  model  as  the  problem;  however  the  failure  on  the  incipient  breaking 
case  which  had  a  steady  non-breaking  solution  remains  unexplained.  For  the  strongly  breaking  waves  of 
Walker  et  al.  (1996),  the  trough  depth  and  the  wavelength  of  the  waves  was  well  predicted  using  RANS, 
but  the  breaking  model  is  too  dissipative  and  leads  to  under-estimation  of  the  following-wave  amplitude. 
For  the  strongest  breaker,  however,  it  was  not  dissipative  enough  and  the  breaker  height  had  to  be 
increased  to  attain  a  non-overturning  solution.  Finally  for  the  cases  of  Furey  et  al.  (2003),  the  best  match 
to  the  data  using  RANS  was  obtained  with  the  free-stream  velocity  set  to  match  the  ADV  velocity 
measurements.  In  general,  the  wavelength  and  the  initial  trough  depth  were  well  predicted. .  The 
breaking  wave  crest  height  was  also  well  predicted  at  the  lower  velocities,  but  over-estimated  at  the 
higher  velocities,  indicating  that  the  level  of  energy  dissipation  caused  by  the  model  is  inappropriate. 

Finally,  the  results  of  two  lower-angle-of-attack  cases  of  Walker  et  al  (1996)  were  used  as  input  to  the 
VSM  radar  backscatter  model  and  radar  cross-section  was  estimated  for  X-Band  HH  and  VV  polarization 
at  45  degrees  incidence,  looking  both  up-wave  and  down-wave.  VSM  assigns  an  empirical  wave 
spectrum  in  the  breaking  region  and  the  calculates  the  propagation  of  the  short  waves  outside  that  region 
using  the  wave-action  balance  equations.  The  main  input  to  VSM  are  the  surface  slopes,  surface 
velocities,  and  a  breaking  mask,  showing  the  extent  of  the  breaking  region.  In  addition  to  RCS,  VSM 
calculates  the  RMS  surface  fluctuation  level,  for  comparison  to  measurements.  The  RMS  surface 
fluctuations  in  all  cases  drop  too  slowly  as  the  waves  propagate  downstream  from  the  breaking  region. 
This  is  traced  to  the  fact  that  the  RANS-predicted  velocity  in  the  breaking  region  is  not  accurate;  this  is  a 
limitation  of  the  breaking  model  which  prescribes  the  surface  velocity  in  the  breaking  region.  Overall  the 
peak  RCS  values  for  the  breaking  region,  tend  to  be  over-estimated  by  up  to  5  dB,  while  the  downstream 
RCS  values  are  over-estimated  by  up  to  25  dB. 

Non-Breaking  Waves 

In  order  to  establish  the  accuracy  of  the  computations  to  be  presented  below  for  breaking  waves,  the  CFD 
codes  were  applied  to  non-breaking  wave  cases.  Computations  were  performed  using  RANS  and 
potential  flow.  Profiles  of  non-breaking  waves  were  measured  by  Duncan  (1983)  and  by  Furey  et  al. 
(2003)  and  are  used  for  comparison.  In  general,  the  results  for  both  approaches  match  the  wavelength  of 
the  experimentally  observed  waves.  The  potential  flow  approach,  due  to  its  lack  of  viscous  effects  will 
always  predict  slightly  larger  lift  on  the  hydrofoil  than  is  realistic  and,  as  a  result,  it  will  over-estimate  the 
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initial  trough  depth  and  the  following  wave  crest  heights.  In  these  results,  the  potential  flow  approach 
predicted  initial  trough  depths  and  crest  height  which  were  generally  10-15%  greater  that  the  RANS 
approach.  For  the  case  of  Furey  et  al.  (2003),  the  potential  flow  approach  over-estimated  the  wave 
amplitude  by  about  10%.  The  RANS  computations  provide  a  better  match  while  slightly  under¬ 
estimating  the  amplitude;  however,  there  is  a  positive  bias  in  the  RANS  results.  For  the  case  of  Duncan 
(1983)  both  approaches  under-estimated  the  wave  amplitude  by  10%  or  more.  The  fact  that  the  measured 
waves  are  larger  than  the  potential  flow  approach  predicts  may  indicate  a  problem  with  the  data. 
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Figure  8.  RANS  computational  results  for  the  flow  over  the  submerged  hydrofoil  moving  with  speed  0.8 
m/s,  depth  of  submergence  0.261m  and  angle  of  attack  5°  (Duncan,  1983):  a)  dynamic 
pressure  field,  b)  axial  velocity  field,  and  c)  turbulence  kinetic  energy  field. 
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The  flow  field  results  for  Duncan’s  case  are  shown  in  Figure  8.  It  can  be  seen  that  the  simulations  predict 
the  expected  flow  behavior.  The  acceleration  of  flow  above  the  hydrofoil  causes  pressure  drop  and  the 
formation  of  the  first  through.  The  surface  then  rises  behind  the  foil  for  the  initial  crest.  In  the  velocity 
field  in  Figure  8b)  and  the  turbulence  kinetic  energy  field  in  Figure  8c)  the  wake  of  the  foil  can  be  seen, 
but  there  is  no  turbulence  near  the  surface  since  the  wave  is  not  breaking.  Figure  9  shows  a  comparison 
of  the  measured  surface-elevation  profile  and  those  computed  using  the  RANS  and  potential  flow 
approaches.  For  both  approaches,  the  wavelength  of  the  waves  in  the  computations  match  the  data  well. 
The  potential  flow  results  under-estimate  the  wave  amplitude  by  about  10%,  and  the  RANS  results  are 
another  10%  lower  that  that.  (The  RANS  results  are  for  a  relatively  low  resolution  computation,  and  as 
discussed  above,  the  waves  behind  initial  crest  exhibit  attenuation;  however,  the  initial  trough  and  crest 
are  unaffected.)  Since  potential  flow  results  would  naturally  over-estimate  the  wave  amplitude,  and  the 
behavior  of  the  RANS  and  potential  flow  results  are  consistent,  there  appears  to  be  a  problem  with  this 
particular  data  set.  This  is  further  supported  by  the  results  shown  above;  this  case  was  used  for  the  grid 
studies,  and  was  thoroughly  examined  using  multiple  grid  resolutions,  topologies,  and  even  different 
RANS  formulations.  All  of  these  were  shown  to  produce  consistent  results. 

Figure  10  shows  a  comparison  of  the  experimental  data  of  Furey  et  al.  (2003)  and  results  from  potential 
flow  and  RANS  computations.  It  should  be  noted  that  the  free  stream  velocity  for  the  calculations  was  set 
to  match  that  measured  using  the  acoustic  anemometer  downstream  of  the  hydrofoil.  For  this  case,  again, 
the  wavelength  of  the  waves  is  well  captured  by  both  approaches.  The  amplitude  of  the  waves  calculated 
using  potential  flow  is  about  10%  larger  than  the  experiments  indicate;  this  is  consistent  with  the 
expectations  for  potential  flow.  The  RANS  results  match  the  amplitude  better,  within  a  few  percent; 
however,  there  is  a  slight  positive  bias  in  the  surface  elevations  estimates. 


Figure  9.  Computational  results  obtained  with  CFDSHIP-IOWA  and  the  potential  flow  code  for  the 
Duncan  0.261  m  case. 
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Figure  10.  Comparison  between  the  results  obtained  with  CFDSHIP-IOWA  and  the  potential  flow  code 
for  the  first  case  by  Furey  et  al.  (2003)  -  CWC  Case  1. 


Breaking  Waves 

Results  for  the  Cases  of  Duncan 

The  measurements  by  Duncan  (1983)  were  performed  using  a  NACA0012  hydrofoil  with  a  chord  of 
0.203  m.  The  upstream  velocity  was  0.8  m/s,  and  the  hydrofoil  angle  of  attack  was  kept  constant  at  5°. 
The  corresponding  Reynolds  and  Froude  numbers  based  on  the  hydrofoil  length  were  1.42xl05  and  0.566. 
Measurements  were  performed  at  varying  submergence  depths,  D.  RANS  computations  were  done  for 
D=0. 1 93  m,  and  D= 0.185  m.  The  first  case  was  called  by  Duncan  ‘incipient  breaking’  because  the  wave 
could  exist  either  as  a  non-breaking  wave  or  a  breaking  wave  depending  on  whether  there  were  outside 
disturbances  to  the  free  surface.  This  case  was  computed  with  and  without  the  wave  breaking  model. 

The  final  case  was  a  stable  breaker,  and  it  was  calculated  using  the  breaking  model. 

Predicted  wave  profiles  together  with  the  measurements  are  shown  in  Figure  1 1 .  The  case  at  D= 0. 1 93  m 
was  run  for  200,000  time  steps  and  did  not  achieve  a  steady  solution;  i.e.  the  wave  train  continued  to 
oscillate  in  amplitude.  After  that  the  simulations  were  restarted  for  another  100,000  time  steps  with  the 
breaking  wave  model  turned  on.  Experience  with  other  breaking  and  non-breaking  cases  indicates  that 
this  should  have  been  sufficient  to  achieve  a  steady-state  solution.  Snap-shots  of  the  wave  field  for  both 
conditions  are  shown  in  Figure  1 1,  where  it  appears  that  the  predicted  amplitude  and  wavelength  are 
reasonable,  but  the  solution  would  not  converge.  Also  shown  in  Figure  1 1  is  the  result  for  the  D=0. 1 85  m 
case.  In  this  case,  similar  behavior  was  observed,  with  no  steady  state  solution  being  obtained.  Again, 
the  wavelength  and  amplitude  predicted  are  reasonable. 
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Figure  11.  Comparison  between  the  predicted  and  measured  wave  profiles  for  the  cases  measured  by 
Duncan  (1983).  Plots  on  the  left  and  bottom  were  obtained  without  the  breaking  wave  model, 
while  plots  on  the  right  were  obtained  with  the  breaking  wave  model  turned  on. 


Flow  field  results  are  shown  in  Figure  12.  for  the  the  case  with  the  hydrofoil  depth  equal  to  0.185  m.  The 
simulation  was  done  with  the  wave  breaking  model.  The  result  of  the  additional  pressure  exerted  by  the 
breaker  on  the  wave  surface  can  be  seen  especially  in  the  velocity  and  turbulence  kinetic  energy  plots. 
Velocity  slows  under  the  breaker.  Figure  12b).  The  velocity  discontinuity  causes  generation  of  turbulence 
that  can  be  seen  in  Figure  12c).  In  these  results,  there  are  no  apparent  problems  which  would  explain  the 
inability  to  achieve  a  steady-state  solution.  Given  that  the  computational  approach  works  well  in  the 
absence  of  wave  breaking,  for  these  breaking  or  incipient  breaking  waves,  the  fault  most  likely  lies  with 
the  wave-breaking  model. 

Results  for  the  Cases  of  Walker  et  al. 

The  measurements  done  by  Walker  et  al.  (1996)  were  performed  using  the  NACA0015  hydrofoil  with  a 
0.304  m  chord.  The  hydrofoil  submergence  depth  was  constant  at  0.267  m.  The  upstream  velocity  varied 
between  1.08  m/s  and  1.10  m/s.  The  corresponding  Reynolds  and  Froude  numbers  were  3.30x10s  and 
0.630.  The  calculations  were  performed  at  three  different  hydrofoil  angles:  0=3°,  0=4°,  and  a=6°.  The 
Walker  et  al.  (1996)  measurements  were  performed  at  considerably  stronger  wave  breaking  conditions 
(more  energy  dissipation)  than  those  of  Duncan  (1983).  These  cases  could  not  be  simulated  without  the 
use  of  the  breaking  wave  model.  Without  the  breaking  model  the  waves  tend  to  overturn  almost 
immediately  after  the  simulation  starts.  Since  the  CFD  model  does  not  admit  overturning  waves,  the 
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simulation  continues  while  the  results  are  not  physical.  Therefore,  all  the  simulations  for  the  Walker  et  al. 
(1996)  cases  were  done  using  the  breaking  wave  model. 


a) 


P 

!  0.0909091 
•0.030303 
•0.151516 
-0,272727 
-0.393939 
-0.515152 

-0.636364 

-0.757576 

•0.878768 


Z 


U 

I  1  17273 
fm  103131 
jM  0  689899 
HI  0  748485 
0  607071 
«  0  465657 
Hi  0  324242 

W  0.162826 

jH  0.0414141 


z 


X 


Figure  12.  Computational  results  for  the  flow  over  the  submerged  hydrofoil  moving  with  speed  0.8  m/s, 
depth  of  submergence  0.185  m  and  angle  of  attack  5°  (Duncan,  1983):  a)  dynamic  pressure 
field,  b)  axial  velocity  field,  and  c)  turbulence  kinetic  energy  field. 


The  results  of  the  RANS  calculations  for  the  cases  of  Walker  et  al.  (1996)  are  shown  in  Figure  13.  Good 
agreement  with  the  experimental  data  can  be  observed  for  the  cases  at  3°  and  4°  hydrofoil  angle.  The 
location  and  depth  of  the  first  trough  is  mostly  in  agreement  with  the  data,  while  the  following  crest  is 
under-predicted.  As  a  result,  the  following  waves  are  also  smaller  than  measured.  The  wavelengths  and 
phase  of  the  waves  are  reasonably  well  predicted.  The  good  agreement  of  the  initial  trough  depths 
indicates  the  basic  RANS  solution  is  reasonable.  The  under-estimation  of  the  breaking  wave  crest  height 
is  a  results  of  the  breaking  wave  model  dissipating  too  much  energy.  It  should  be  noted  that  the  breaking 
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model  was  mainly  developed  using  the  data  of  Duncan  (1983)  for  waves  which  were  not  nearly  as 
dissipative  as  these. 


h'=0.640H 


h*=0.640H 


Figure  13.  Comparison  between  the  predicted  and  measured  wave  profiles  for  the  cases  measured  by 
Walker  et  al.  (1996).  Simulations  were  performed  at  different  breaker  heights,  h*,  used  in  the 
breaking  wave  model. 


For  the  6°  case,  shown  at  the  bottom  of  Figure  13,  the  simulation  became  unstable  and  the  wave  started  to 
overturn,  leading  to  a  non-physical  solution.  The  reason  for  this  was  that  the  breaking  wave  model  was 
not  able  to  dissipate  enough  energy  to  stabilize  the  wave.  For  this  case,  the  experimentally  observed 
breaker  height  h*  was  approximately  equal  to  the  wave  height  H.  Setting  h*=H  (rather  than  the  literature 
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value  of  0.64//)  resulted  in  a  stable  calculation,  but  the  breaking  region  was  now  confined  to  a  small 
region  in  the  trough  of  the  wave.  This  results  in  a  strange  inflection  in  the  trough  of  the  wave  and  is 
probably  non-physical. 

Results  for  the  Cases  of  Furey  et  al. 

Recently,  Furey  et  al.  (2003)  performed  wave  measurements  in  considerably  larger  channel  than  the  two 
previous  experiments.  They  used  hydrofoil  NACA0012  profile  with  a  chord  1.8  m  in  length.  The 
distance  between  the  channel  floor  and  the  hydrofoil  attachment  point  was  0.9  m.  Submergence  depth 
varied  slightly  between  1.75  m  and  1.79  m.  Hydrofoil  angle  of  attack  was  set  to  5.5+0.50.  Measurements 
were  performed  at  different  inlet  velocities  between  1 .8  m/s  and  2.6  m/s  that  were  also  used  in  the 
calculations:  w=1.8  m/s  and  £ts=5.5°,  u-2.2  m/s  and  ot=5.50',  «=2.3  m/s  and  at=5.5°;  w=2.6  m/s  and 
<*=5.5°. 

These  inlet  velocities  were  measured  with  a  pitot  tube  upstream  of  the  foil.  The  simulations  were 
performed  using  the  breaking  wave  model  with  the  original  parameters  (/i*=0.64/7).  Comparison  between 
the  predictions  and  the  measurements  is  shown  in  the  left  column  of  Figure  14.  The  agreement  is  the 
worst  for  the  first,  non-breaking  case.  The  agreement  slightly  improves  at  higher  inlet  velocities  and 
under  breaking  wave  conditions.  However,  wave  amplitudes  are  under  predicted  for  all  the  cases.  All 
calculated  troughs  are  located  upstream  of  their  measured  location. 

Surface  velocity  was  measured  also  with  the  acoustic  doppler  velocimeter  (ADV).  The  ADV-  measured 
velocities  were  higher  than  the  values  obtained  using  the  pitot  tube.  Another  set  of  calculations  was 
performed  with  these  values  at  the  inlet  and  the  results  are  shown  in  the  right  column  of  Figure  14.  The 
agreement  between  the  data  and  the  measurements  are  considerably  better.  The  predicted  location  of 
troughs  and  crests  is  in  agreement  with  the  measurements.  The  overall  wave  profiles  are  in  excellent 
agreement  for  the  first  two  cases.  Further  increase  in  the  inlet  velocity  causes  slightly  under-predicted 
troughs  and  over  predicted  crests.  A  possible  cause  for  such  behavior  might  be  the  deficiency  of  the  wave 
breaking  model. 
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Figure  14.  Comparison  between  the  predicted  and  measured  wave  profiles  for  the  data  by  Furey  et  al 
(2003).  Simulations  were  performed  at  inlet  velocities  measured  by  the  Pitot  tube  (left 
column)  and  by  the  Acoustic  Doppler  Velocimeter  (right  column). 
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Radar  Backscatter  Modeling 

For  the  cases  of  Walker  et  al.  (1996)  detailed  radar  backscatter  measurements  are  available,  as  presented 
in  that  study,  and  later  in  Ericson  et  al.  (1999).  The  modeling  approach  to  be  used  is  the  small- 
perturbation  method,  or  SPM,  advanced  by  Valenzuela  (1978).  This  method  was  applied  to  the  situation 
examined  here  by  Ericson  et  al.  (1999),  using  the  experimentally  measured  short-wave  spectrum  and 
surface  slope  as  input  to  the  model  The  results  are  shown  in  Figure  15,  where  it  is  seen  that  the 
agreement  is  quite  good  for  the  region  downstream  of  the  breaking  wave  crest.  The  backscatter  from  the 
breaking  crest  is  not  well  predicted  using  this  model,  since  the  surface  roughness  exceeds  the  limits  of 
applicability  of  the  theory.  For  the  up-wave  look  directions  the  breaking  crest  returns  are  over-predicted 
for  both  HH  and  VV,  while  for  down-wave  look  direction  the  breaking  crest  return  is  slightly  over¬ 
predicted  for  VV  and  under-predicted  for  HH.  In  general,  the  lack  of  polarization  dependence  of  the 
scattering  in  the  breaking  region  is  not  captured. 


Position  Along  Wavetank  (m)  Position  Along  Wavetank  (m) 

(c)  (d) 

Figure  15.  Radar  cross-section  results  for  the  small-perturbation  method  (SPM)  calculated  using  the 
experimentally  measured  short-wave  spectrum,  compared  to  experimental  observations  for 
both  HH  and  VV  polarization  (from  Ericson  et  al.  1999):  a)  3  degree  angle  of  attack,  down- 
wave  look  direction;  b)  3  degree  angle  of  attack,  up-wave  look  direction;  c)  4  degree  angle  of 
attack,  down-wave  look  direction;  d)  4  degree  angle  of  attack,  up-wave  look  direction. 


23 


DW200392.8.F.Nov04 


It  should  be  noted  that  in  Ericson  et  al.  (1999),  the  integral  equation  method  (IEM)  was  the  best¬ 
performing  model  over  the  entire  region  of  interest.  For  this  study,  the  VSM  code  was  used  for  the 
backscatter  calculations.  The  IEM  implementation  in  VSM  is  only  approximate,  and  not  the  complete 
model  advanced  by  Fung  et  al.  (1992)  and  used  by  Ericson  et  al.  (1999).  For  this  reason  the  comparisons 
presented  below  are  based  on  the  SPM,  rather  than  IEM.  Even  though  the  agreement  of  SPM  model  is 
not  exact,  valid  comparisons  can  still  be  made.  If  the  hydrodynamic  modeling  produces  a  mean  surface 
elevation  profile  and  a  short  wave  spectrum  which  matches  the  experimental  spectrum  observations,  the 
RCS  predictions  will  match  those  shown  in  Figure  15. 

The  hydrodynamic  results  generated  using  the  RANS  models  shown  above  were  used  as  input  for  VSM. 
The  RCS  results  from  VSM  for  the  3  degree  angle  of  attack  are  shown  in  Figure  16  along  with  the 
predicted  mean  surface  elevation  and  RMS  surface  fluctuation  level.  Figure  16  also  includes 
experimental  observations  for  all  these  quantities. 

In  Figure  16a)  the  predicted  mean  surface  elevation  is  compared  to  data.  As  discussed  above,  the  height 
of  the  breaking  wave  crest,  and  that  of  the  following  crests  was  under-estimated.  The  model  predictions 
for  the  RMS  surface  elevation  were  obtained  by  integrating  the  short-wave  spectrum  estimated  using 
VSM.  Two  results  are  shown  in  Figure  16b)  for  comparison  to  the  experimental  data:  One  is  that 
obtained  when  VSM  is  allowed  to  determine  the  breaking  region  based  on  the  vertical  acceleration.  The 
other  is  for  the  breaking  region  defined  by  the  RANS  wave-breaking  model.  For  the  former  case,  the 
peak  RMS  surface  elevation,  which  occurs  in  the  breaking  crest  region,  is  under-estimated  by  roughly  a 
factor  of  two.  The  modeled  RMS  decays  more  slowly  than  the  experiments  indicate,  and  the  agreement  is 
not  bad  downstream  of  the  breaking  region.  Since  the  short  waves  which  establish  the  RMS  in  this 
downstream  region  are  generated  in  the  breaking  region,  and  their  energy  level  is  too  low  in  the  breaking 
region,  the  agreement  downstream  is  clearly  the  result  of  offsetting  errors  in  the  generation  and 
propagation  of  the  short  waves.  When  the  RANS  code  is  used  to  determine  the  breaking  region,  the  RMS 
is  roughly  twice  the  experimental  value  in  the  breaking  region,  and  again  decays  too  slowly  downstream. 
For  this  case,  the  RMS  remains  elevated  over  the  following  waves.  The  reason  for  the  under-estimation 
for  the  decay  in  the  RMS  surface  elevation  is  discussed  below. 

Figure  16c)-  Figure  16f)  show  VSM  predictions  of  RCS  for  the  3  degree  angle  of  attack  case  for  HH  and 
VV  polarization  and  for  up-wave  and  down-wave  look  direction.  Two  computational  results  are  shown, 
one  using  VSM  to  determine  the  breaking  region  and  the  other  using  the  RANS-indicated  breaking 
region,  along  with  experimental  data  from  Ericson  et  al.  (1999).  First  we  examine  the  results  for  the 
RANS-indicated  breaking.  The  computational  radar  backscatter  results  for  the  breaking  wave  crest  for 
HH  up-wave  (Figure  16c)  match  the  experimental  data,  while  for  HH  down-wave  (Figure  16d)  the 
predictions  are  roughly  5  dB  high.  Figure  15a)  shows  that  the  ‘correct’  SPM -predicted  down-wave  HH 
backscatter  from  the  breaking  crest  is  actually  lower  by  more  than  10  dB  than  for  the  up-wave  look 
direction.  Figure  16c)  and  d)  shows  that  the  RANS-based  up-wave  and  down-wave  returns  are  nearly 
equal.  This  discrepancy  results  from  both  the  error  in  the  size  of  the  breaking  region  and  the  lower- 
amplitude,  relatively  flat  breaking  crest.  The  computational  results  for  VV  up-wave  in  the  breaking  crest 
region  (Figure  16e)  are  a  few  dB  higher  than  the  HH  results,  which  is  consistent  with  the  data  and  with 
the  behavior  of  the  SPM  model.  Similar  behavior  is  seen  in  the  W  down-wave  predictions  for  the 
breaking  crest  (Figure  16f).  Downstream  of  the  breaking  crest,  the  decay  of  the  radar  cross-section 
mirrors  the  too-slow  decay  of  the  RMS  surface  roughness,  yielding  over-estimations  of  the  backscatter  by 
as  much  as  25  dB,  while  the  modulations  have  about  the  right  amplitude.  For  the  VSM-indicated 
breaking,  the  results  are  generally  similar;  however,  the  lower  RMS  fluctuation  levels  results  in  generally 
lower  backscatter. 
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Figure  16.  Hydrodynamic  and  SPM  radar  backscatter  predictions  for  3  degree  case  of  Walker  et  al. 

(1996):  a)  mean  surface  elevation;  b)  RMS  surface  elevation;  c)  radar  cross  section  for  HH 
polarization,  up-wave  look  direction;  d)  HH,  down-wave;  e)  VV,  up-wave;  f)  VV,  down- 
wave 
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The  lower  predicted  decay  rate  of  the  short  waves  downstream  of  the  breaking  region  shown  in  Figure 
16b)  is  related  to  the  prediction  of  the  velocity  field  by  the  RANS  models.  As  the  waves  generated  in  the 
breaking  region  propagate  downstream  they  interact  with  an  accelerating  fluid  velocity.  This  acceleration 
will  reduce  the  wave  energy,  i.e.  the  waves  are  ‘stretched’  as  they  leave  the  breaking  region,  lowering 
their  amplitude.  The  lack  of  a  rapid  drop  in  RMS  just  downstream  of  the  breaking  region  indicates  that 
the  waves  do  not  undergo  this  acceleration  in  the  model  calculations.  This  can  come  about  in  two 
possible  ways.  Since  the  wave  spectrum  is  prescribed  in  the  indicated  breaking  region,  it  is  is  only 
influenced  by  the  fluid  velocities  outside  this  region;  hence,  if  the  indicated  breaking  region  extends 
beyond  this  area  of  rapid  acceleration,  the  wave  will  never  experience  it.  Alternatively  if  the 
hydrodynamic  model  for  breaking  under-estimates  this  acceleration,  the  waves  will  not  be  reduced  in 
amplitude  as  much  either.  The  latter  explanation  fits  the  present  case;  Figure  17  shows  a  comparison  of 
experimental  data  for  the  surface  velocity  from  Ericson  et  al.  (1999)  to  that  obtained  using  RANS.  In  the 
experiments,  the  velocity  in  the  breaking  region  is  essentially  zero;  however,  the  RANS  calculations 
indicate  that  the  velocity  in  the  breaking  region  is  only  reduced  by  half.  This  indicates  that  the  breaking 
model  is  producing  a  poor  estimate  of  the  velocities  in  the  breaking  region  and  this  directly  affects  the 
short  wave  spectrum  prediction  and  that  for  the  RMS  surface  elevation. 


0.5  1  1.5  2  2.5 

x(m) 

Figure  17.  Comparison  of  measured  (upper)  and  predicted  (lower)  streamwise  velocity  components  at 
surface  for  the  3  degree  angle  of  attack  wave  from  Ericson  et  al.  (1999).  Experimental  data 
from  Ericson  et  al.,  computations  using  CFDSHIP-IOWA.  The  dashed  red  line  indicates  the 
location  of  breaking  in  the  RANS  predictions. 
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Figure  18.  Hydrodynamic  and  SPM  radar  backscatter  predictions  for  4  degree  case  of  Walker  et  al. 

(1996):  a)  mean  surface  elevation;  b)  RMS  surface  elevation;  c)  radar  cross  section  for  HH 
polarization,  up-wave  look  direction;  d)  HH,  down-wave;  e)  VV,  up-wave;  f)  VV,  down- 
wave. 
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Figure  18  shows  the  results  from  VSM  for  the  4  degree  case  of  Walker  et  al.  Overall  the  behavior  is 
similar  to  that  of  the  3  degree  case  shown  in  Figure  16.  The  main  difference  is  the  small  degree  of 
modulation  in  the  downstream  RCS  for  RANS-predicted  breaking  region.  This  is  a  result  of  the  lower 
amplitude  of  the  following  waves,  when  compared  to  the  3  degree  case,  consistent  with  greater  energy 
dissipation  due  to  breaking  than  in  the  3  degree  case.  The  difference  between  the  downstream  RCS 
modulations  for  the  RANS-  and  VSM-indicated  breaking  regions  is  due  to  the  broader  spectrum  (from  the 
larger  breaking  region)  in  the  RANS  case,  and  the  smaller  influence  of  the  orbital  velocities  on  the  longer 
waves. 
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Conclusions  and  Recommendations 

Conclusions 

This  study  had  two  goals,  the  first  of  which  was  to  examine  approaches  for  calculating  the  hydrodynamics 
of  breaking  waves  and  to  assess  the  impact  of  hydrodynamic  model  errors  on  the  prediction  of  radar 
backscatter.  To  address  this,  RANS  computations  of  stationary  hydrofoil-generated  breaking  waves  were 
carried  out,  including  the  modeling  of  the  breaking  region.  These  results  were  compared  to  experimental 
data.  A  subset  of  these  hydrodynamic  results  was  used  as  input  to  the  SPM  radar  backscatter  model  in  the 
VSM  code  and  the  RCS  predictions  were  compared  to  those  obtained  using  the  same  model  with 
experimentally  measured  hydrodynamic  inputs.  The  second  goal  of  this  study  was  to  use  the  information 
gained  through  this  effort  to  define  the  research  needs  in  this  area. 

The  results  show  that  in  order  to  obtain  a  RANS  solution,  the  domain  must  be  large  and  the  grid 
sufficiently  stretched  to  eliminate  reflection  from  the  upstream  and  downstream  domain  boundaries. 
Additional  damping  of  the  waves  is  sometimes  needed.  What  was  also  learned  is  that  as  long  as  the  flow 
near  the  foil  is  well  resolved,  the  near-field  solution  for  the  free  surface  elevation  and  velocity  field  will 
be  accurately  captured;  however,  if  the  grid  resolution  is  not  adequate,  the  far  field  waves  will  be 
attenuated.  This  near-field  accuracy  vs.  far-field  attenuation  has  been  seen  in  Kelvin-wake  predictions, 
and  is  tied  directly  to  grid  resolution  and  its  impact  on  wave  propagation. 

For  the  non-breaking  waves,  the  main  results  indicate  that  while  RANS  and  potential  flow  approaches  are 
capable  of  accurately  calculating  the  wavelength  of  hydrofoil-generated  stationary  breakers,  potential 
flow  methods  tend  to  over-estimate  the  wave  amplitude,  while  RANS  appears  to  under-estimate  the 
amplitude.  It  should  be  noted  that  the  amount  of  data  available  for  comparison  is  limited  in  both  extent, 
and  type.  Velocity  fields  along  with  surface  elevation  profiles  would  allow  more  detailed  comparisons 
and  more  complete  validation  of  calculations. 

For  the  breaking  wave  cases,  the  breaker  model  was  shown  to  have  qualitatively  correct  behavior;  i.e. 
energy  is  dissipated  and  turbulence  is  produced  in  the  breaking  region  and  the  wave  amplitude  is  reduced, 
but  was  not  robust  or  particularly  accurate.  There  were  cases  which  could  not  be  driven  to  a  converged 
solution.  For  the  strongly  breaking  waves  of  Walker  et  al.  (1996),  the  trough  depth  and  the  wavelength  of 
the  waves  was  well  predicted  using  RANS,  but  the  breaking  model  was  too  dissipative  and  led  to  under¬ 
estimation  of  the  following-wave  amplitude.  For  the  strongest  breaker,  however,  it  was  not  dissipative 
enough  and  the  breaker  height  had  to  be  increased  to  attain  a  non-overturning  solution.  For  the  cases  of 
Furey  et  al.  (2003),  the  wavelength  and  the  initial  trough  depth  were  well  predicted. .  The  breaking  wave 
crest  height  was  also  well  predicted  at  the  lower  velocities,  but  over-estimated  at  the  higher  velocities, 
indicating  that  the  level  of  energy  dissipation  caused  by  the  model  is  inappropriate  at  higher  speeds. 

Finally,  the  results  of  two  lower-angle-of-attack  cases  of  Walker  et  al  (1996)  were  used  as  input  to  the 
VSM  radar  backscatter  model  and  radar  cross-section  was  estimated  for  X-Band  HH  and  VV  polarization 
at  45  degrees  incidence,  looking  both  up-wave  and  down-wave.  In  addition  to  RCS,  VSM  calculates  the 
RMS  surface  fluctuation  level,  for  comparison  to  measurements.  The  RMS  surface  fluctuations  in  all 
cases  drop  too  slowly  as  the  short  waves  propagate  downstream  from  the  breaking  region.  This  is  traced 
to  the  fact  that  the  RANS-predicted  velocity  in  the  breaking  region  is  not  accurate;  this  is  a  limitation  of 
the  breaking  model  which  prescribes  the  surface  velocity  in  the  breaking  region.  Overall  the  peak  RCS 
values,  for  the  breaking  region,  tend  to  be  over-estimated  by  as  much  as  12  dB,  and  the  downstream  RCS 
values  tend  to  be  over-estimated  by  as  much  as  25  dB.  Hence,  hydrodynamic  errors  incurred  in  modeling 
breaking  waves  can  lead  to  significant  errors  in  radar  backscatter  estimates. 
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Research  Needs 


In  order  to  predict  the  behavior  of  breaking  waves  for  flows  of  practical  interest,  approaches  such  as  the 
Reynolds-averaged  Navier-Stokes  equations  must  be  employed.  This  type  of  approach  requires  that 
some  of  the  details  of  the  flow  must  be  represented  by  models  which  need  to  accurately  reflect  the 
physical  mechanisms  at  play.  The  results  of  this  study  indicate  areas  where  more  knowledge  and 
understating  are  required  in  order  to  develop  the  needed  models,  and  point  to  a  requirement  for  new,  more 
detailed  information  on  the  behavior  of  breaking  waves.  There  are  two  ways  in  which  this  information 
can  be  obtained:  The  first  is  through  direct  numerical  simulations  (DNS)  based  on  numerical  solution  of 
the  exact  governing  equations  for  the  three-dimensional  time-evolution  of  breaking  waves.  DNS  allows 
the  details  of  the  waves  to  be  examined  both  qualitatively  and  quantitatively,  and  can  provide  access  to 
un-measurable  quantities  such  as  the  subsurface  pressure  fluctuations.  These  detailed  data  can  be  used  to 
guide  the  development  of  the  needed  models.  The  drawback  to  this  approach  is  that  presently  it  can  only 
be  applied  to  small-scale  waves  (i.e.  waves  substantially  less  than  a  meter  in  wavelength)  and  even  then  it 
taxes  available  computing  resources.  The  second  approach  is  experimental  studies  of  breaking  waves, 
where  a  more  limited  set  of  information  can  be  obtained  using  available  measurement  technology,  but  this 
can  be  done  at  larger  scale,  approaching  the  scale  of  situations  of  practical  interest.  Detailed 
measurement  of  surface  structure,  velocity  and  subsurface  turbulence  are  needed.  Through  a  coordinated, 
combined  approach  of  DNS  and  experiments,  the  information  necessary  to  develop  first-principles 
models  for  wave  breaking  and  the  associated  small-scale  disturbances  important  for  radar  backscatter  can 
be  obtained. 

Beyond  the  general  need  for  more  detailed  information  about  the  behavior  of  breaking  waves,  the  results 
of  this  study  indicate  some  specific  areas  where  improvements  must  be  made.  The  first  is  application  of 
computational  models  to  steep  waves,  where  more  robust  and  efficient  computational  approaches  are 
clearly  needed.  The  second  is  modeling  of  the  effects  of  wave  breaking,  where  issues  such  as  modeling 
the  turbulence  in  the  breaking  region,  and  appropriate  ways  to  represent  free-surface  boundary  conditions, 
etc.  in  the  breaking  region  need  to  be  addressed.  And  the  third  is  use  of  hydrodynamic  model  results  in 
radar  backscatter  prediction.  These  point  clearly  to  research  needs  in  these  areas  and  are  presented  below. 

A  final  issue  is  the  framework  within  which  modeling  of  breaking  waves  should  be  addressed.  Unlike 
potential-flow  methods,  Reynolds-averaged  approaches  do  not  exclude  any  of  the  physical  processes 
which  naturally  occur  in  wave  breaking;  however,  potential-flow  methods  are  significantly  less 
demanding  computationally.  While  RANS  approaches  should  be  able  to  predict  wave  breaking 
accurately,  given  appropriate  models,  it  would  be  desirable  to  be  able  to  use  potential  flow  methods.  This 
is  also  discussed  below. 

Application  of  CFD  to  Steep  Waves 

One  of  the  unexpected  developments  in  this  study  was  the  difficulty  encountered  in  carrying  out  the 
computations  of  apparently  simple,  hydrofoil-generated  waves,  even  for  steep,  non-breaking  cases.  One 
of  the  main  difficulties  was  due  to  the  lack  of  appropriate  radiation  boundary  conditions,  necessary 
eliminate  wave  reflections  from  the  domain  boundaries.  Present  approaches  extend  the  size  of  the  domain 
and  introduce  either  implicit  or  explicit  damping  of  the  waves,  so  that  the  waves  never  actually  reach  the 
boundaries.  This  leads  to  unreasonably  large  computational  domains  and  significantly  greater 
computational-time  requirements.  The  other  main  issue  was  the  need  for  robustness  in  free-surface 
computations.  The  free-surface-tracking  approaches  implemented  in  the  codes  used  above,  where  the 
computational  domain  is  mapped  to  the  instantaneous  free  surface,  implicitly  assumes  that  the  surface  is  a 
single- valued  function  of  the  horizontal  coordinates.  This  eliminates  the  possibility  of  overturning  in  the 
computations  and  so  if  breaking  is  initiated  by  overturning,  this  cannot  be  captured  directly.  A 
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computational  approach  which  admits  a  multiple-valued  free  surface  description,  such  as  a  level-set  or 
more  robust  surface-tracking  approach  is  needed. 

Modeling  the  Effects  of  Wave  Breaking 

If  a  wave  becomes  steep  enough,  a  portion  of  its  energy  is  dissipated  via  wave  breaking.  Overturning  of 
the  crest  and  turbulence  generation,  accompanied  by  transfer  of  momentum  from  the  wave  orbital 
motions  to  the  mean  drift  flow  accomplish  this  dissipation.  The  breaker  model  investigated  here 
approaches  wave  breaking  by  attempting  to  add  an  appropriate  amount  of  pressure  to  the  free  surface  in 
the  breaking  region,  in  an  attempt  to  stabilize  it  and  keep  it  from  overturning.  The  rate  of  work  done  on 
the  surface  by  this  pressure  equals  the  dissipation  due  to  breaking.  This  approach  ignores  the  fact  that  in 
even  in  the  transient  leading  to  ultimately  steady  breaking,  some  overturning  occurs,  and  it  requires 
significant  empiricism  in  specifying  that  stabilizing  pressure.  A  final,  robust  model  for  wave  breaking 
will  probably  not  be  of  this  form.  The  results  of  this  study  raises  additional  significant  issues  which  will 
need  to  be  addressed  by  the  ultimate  approach.  These  include  1)  determining  the  onset  of  wave  breaking, 
2)  modeling  the  generation  of  surface  disturbances  in  the  breaking  region,  and  3)  modeling  the  generation 
of  waves  in  the  breaking  region  and  their  subsequent  propagation. 

Modeling  approaches  for  the  onset  of  wave  breaking  do  not  appear  to  replicate  observed  behavior  very 
well.  As  a  wave  becomes  steeper,  the  water  particle  velocities  increase  roughly  in  proportion  to  its 
steepness.  Eventually,  as  the  steepness  increases,  the  water  particles  near  the  crest  begin  to  move  at  a 
speed  greater  that  the  wave  propagation  velocity.  This  leads  to  overturning  and  collapse  of  the  wave  crest 
into  breaking.  In  present  modeling  approaches,  the  onset  of  breaking  is  identified  by  examining  wave 
steepness  (greater  than  15  degrees)  or  downward  accelerations  (greater  than  0.1 5g).  The  problem  with 
these  approaches  is  that  the  thresholds  are  generally  required  to  be  substantially  lower  than  observed  in 
waves  just  before  breaking  occurs  naturally.  In  the  final  model  for  wave  breaking,  the  natural  evolution 
of  a  steep  wave  to  a  state  of  breaking  should  be  captured. 

The  generation  of  surface  disturbances  in  the  breaking  region  is  important  for  radar  backscatter  as  well  as 
for  hydrodynamics.  The  effects  on  radar  backscatter  are  clear.  For  hydrodynamics,  the  disturbances  are 
caused  by  the  subsurface  turbulent  fluctuations,  and  this  process  affects  the  dynamics  of  the  turbulence 
and  thereby  the  overall  evolution  and  behavior  of  the  flow.  The  generation  of  surface  disturbances  has 
not  been  investigated  significantly,  and  it  is  not  clear  how  modeling  of  this  can  be  accomplished  in  a 
RANS  context.  Presently  these  disturbances  are  treated  as  waves  with  a  prescribed  spectrum,  but  it  is  not 
clear  that  all  the  surface  fluctuations  in  the  breaking  region  can  be  appropriately  described  this  way.  This 
is  an  area  where  guidance  from  DNS  computations  of  breaking  waves  could  be  of  great  use  for  model 
development. 

The  disturbances  described  above,  while  important  in  their  own  right,  also  serve  as  initial  conditions  for 
waves  which  propagate  through  and  beyond  the  breaking  region  and  over  the  water  surface.  These  waves 
contribute  significantly  to  radar  backscatter.  They  also  interact  with  the  surrounding  velocity  field  and 
can  steepen  and  break  themselves.  These  waves  are  described  statistically  by  the  wave  action  balance 
equations  but,  given  the  strong  interactions  of  the  waves  with  other  disturbances  in  the  breaking  region 
and  the  large  velocity  gradients  involved,  it  is  not  clear  that  this  description  which  assumes  slow  spatial 
and  temporal  changes,  is  adequate.  In  addition,  there  are  source  terms  (turbulent  wave  generation  and 
dissipation)  which  are  recognized  as  necessary  which  are  not  included  in  the  present  formulations  and 
require  development.  This  is  again  an  area  where  DNS  would  be  instructive  and  useful. 
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Hydrodynamic  Results  for  Radar  Backscatter  Modeling 

Modeling  radar  backscatter  presently  relies  primarily  on  a  statistical  description  of  the  small-scale 
disturbances  on  the  water  surface  along  with  information  on  the  mean  surface  slope.  In  an  ideal  world,  a 
three-dimensional  time-varying  description  of  the  water  surface  at  all  scales  would  be  obtainable,  and 
from  this  a  radar  backscatter  model  could  make  use  of  whatever  statistical  description  of  the  surface  is 
most  appropriate.  Since  this  detailed  description  is  not  obtainable  for  practical  flows,  hydrodynamic 
models  which  produce  the  necessary  statistical  description  should  be  the  focus  of  future  research.  As  our 
understanding  of  the  behavior  of  the  water  surface  in  the  breaking  region  becomes  more  clear,  it  may 
prove  that  the  wave  spectrum  alone  is  an  inadequate  description.  In  that  case,  alternative  descriptions 
would  need  to  be  explored  and  may  influence  the  modeling  approaches  used  for  the  small-scale 
disturbances,  and  possibly  even  those  adopted  for  radar  backscatter  modeling. 

Computational  Framework  for  Wave  Breaking  Models 

A  final  issue  is  the  framework  within  which  modeling  of  breaking  waves  should  be  addressed.  The  exact 
equations  are  known,  but  cannot  be  solved  for  large-scale  flows  and  so  reduced  forms  of  the  equations 
must  be  used.  A  Reynolds  averaged  approach,  while  approximate,  makes  no  assumptions  about  the 
nature  of  the  flow,  and  includes  turbulence  and  viscous  effects.  Potential  flow  formulations  require  that 
the  flow  be  inviscid  and  irrotational,  and  exclude  turbulence.  Since  turbulence  and  viscous  effects  appear 
to  play  significant  roles  in  breaking  waves,  the  RANS  approach  would  appear  to  be  a  natural  framework 
in  which  to  treat  breaking  waves;  however,  the  speed  at  which  potential  flow  computations  can  be  carried 
out  make  this  an  attractive  platform  for  breaking  wave  computations.  Their  usefulness  will  hinge  on 
development  of  a  physically  reasonable  approach  for  modeling  wave  breaking.  While  potential  flow 
approaches  may  ultimately  prove  to  be  too  empirical  or  ad-hoc,  they  warrant  investigation. 


32 


DW200392.8.F.Nov04 


References 


von  Brummelen,  E.H.,  Raven,  H.C.,  Koren,  B.,  “Efficient  numerical  solution  of  steady  free-surface 
Navier-Stokes  flows, ,”  J.  Comp.  Physics,  174,  pp.  120-137,  2001. 

Cointe,  R.,  Tulin,  M.P.,  “A  theory  of  steady  breakers,”  J.  Fluid  Mechanics,  276,  pp.  1-20,  1994. 

Duncan,  J.H.,  ‘The  breaking  and  non-breaking  wave  resistance  of  a  two-dimensional  hydrofoil,”  J.  Fluid 
Mechanics,  126,  pp.  507-520,  1983. 

Ericson,  E.A.,  Wackerman,  C.C.  The  theory  and  implementation  of  ESM,  ERIM  International,  Inc. 
Report  No.  1001 125-1-T,  2000. 

Ericson,  E.A.,  Lyzenga,  D.R.,  Walker,  D.T.  “Radar  Backscatter  from  Stationary  Breaking  Waves,”  J. 
Geophys.  Res.,  104,  pp.  29,679-29,695,  1999. 

Fung,  A.  K.,  Z.  Li,  and  K.  S.  Chen,  Backscattering  from  a  randomly  rough  dielectric  surface,  IEEE  Trans. 
Geosci.  Remote  Sens.,  30,  356-369,  1992. 

Furey,  D.A.,  Fu,  T.C.,  Karion,  A.,  Sur,  T.W.,  Rice,  J.R.,  Walker,  D.C.,  “Experimental  Study  of  the  Wave 
Field  Produced  by  Submerged  Hydrofoil,”  8th  Int.  Conf.  on  Numerical  Ship  Hydrodynamics,  Busan, 
Korea,  September,  2003. 

Iafrati,  A.,  Di  Mascio,  A.,  Campana,  E.F.,  “A  level  set  technique  applied  to  unsteady  free  surface  flows,” 
Int.  J.  Num.l  Methods  Fluids,  35,  pp.  281-297, 2001. 

Komen,  G.J.,  Cavaleri,  L.,  Donelan,  M.,  Hasselmann,  K.,  Hasselmann,  S.,  Janssen,  P.A.E.M.,  Dynamics 
and  Modelling  of  Ocean  Waves,  Cambridge,  1994. 

Kuethe,  A.  M.,  Chow,  C.,  Foundations  of  Aerodynamics:  Bases  of  Aerodynamic  Design,  4th  ed.  John 
Wiley  &  Sons,  Inc.,  1986. 

Muscari,  R.,  Di  Mascio,  A.,  “A  model  for  the  simulation  of  steady  spilling  breaking  waves,”  J.  Ship 
Research,  47,  pp.  13-23,  2003. 

Paterson,  E.G.,  Wilson,  R.V.,  Stem,  F.,  “General-Purpose  parallel  Unsteady  RANS  Ship  Hydrodynamics 
Code:  CFDSHIP-IOWA,”  HHR  Report  No.  432,  The  University  of  Iowa,  November  2003. 

Rhee,  S.H.,  Stem,  F.,  “RANS  Model  for  Spilling  Breaking  Waves,”  J.  Fluids  Engineering,  124,  pp.  424- 
432,  2002. 

Scullen,  D.  C.,  Tuck,  E.  O.,  “Nonlinear  free-surface  flow  computations  for  submerged  cylinders,” 

Journal  of  Ship  Research  39, 185-193,  1995. 

Scullen,  D.  C.,  “Accurate  computation  of  steady  nonlinear  free-surface  flows,”  Ph.D.  Thesis,  The 
Department  of  Applied  Mathematics,  The  University  of  Adelaide,  1998. 

Valenzuela,  G.R.,  “Theories  for  the  interaction  of  electromagnetic  and  oceanic  waves:  A  review,” 
Boundary  Layer  Meteorol.  13,61-85,  1978. 


33 


DW200392.8.F.Nov04 


Walker,  D.T.,  Lyzenga,  D.R.,  Ericson,  E.A.,  Lund,  D.E.,  “Radar  Backscatter  and  Surface  Roughness 
Measurements  for  Stationary  Breaking  Waves,”  Proc.  R.  Soc.  Lond.,  452,  pp.  1953-1984,  1996. 


34 


DW200392.8.F.Nov04 


Appendix:  Nonlinear  Free-Surface  Potential-Flow  Computation 

Involving  a  Hydrofoil 


Methodology 

The  two-dimensional  potential  flow  and  nonlinear  free-surface  that  results  from  submerging  a  hydrofoil 
some  fixed  depth  beneath  the  free  surface  of  a  finite  depth  of  water  was  computed.  The  developed 
computational  method  is  a  combination  of  two  potential-flow  approaches.  An  iterative  de-singularized 
approach  [1,  2]  is  used  to  compute  the  nonlinear  free-surface,  and  a  linearly- varying  vortex-strength  panel 
approach,  similar  to  common  vortex-panel  methods  [3, 4],  is  used  to  compute  the  flow  about  the 
hydrofoil.  Combining  these  two  approaches  was  achieved  through  three  steps. 

First,  the  de-singularized  nonlinear  free-surface  potential-flow  approach,  presented  by  Scullen  [1,2],  was 
implemented  and  verified.  This  approach  utilizes  a  distribution  of  discrete  singularities  (i.e.  sources, 
doublets,  and  vortices),  placed  outside  the  fluid  domain,  and  an  iterative  approach  to  simultaneously 
satisfy  the  various  boundary  conditions  at  the  free  surface.  At  each  iteration,  new  singularity  strengths  are 
determined  by  the  solution  of  a  system  of  linear  equations,  and  the  free-surface  elevation  is  updated  based 
on  the  new  pressure  distribution  at  the  free  surface. 

The  following  set  of  equations  describes  the  flow  that  is  considered: 


V20  =  o  in  the  fluid  domain, 

(1) 

E  =  — ^  =  0  and 

(2) 

Dt 

p  =  0  on  the  free  surface, 

(3) 

V  (/>  — »  Um  upstream,  and 

(4) 

V<p  •  n  =  0  on  the  body. 

(5) 

Here,  $  represents  the  velocity  potential,  p  is  the  pressure  in  excess  of  atmospheric,  C/„  is  the  upstream 
speed  of  the  fluid,  and  n  is  the  unit-normal  vector  pointing  out  of  the  body.  Since  the  pressure  is 
analytically  determined,  via  Bernoulli’s  equation,  as  a  function  of  the  fluid  depth  and  velocity  potential 
and  (/)  identically  satisfies  Equation  1,  the  remaining  task  is  to  compute  the  velocity  potential  and  free- 
surface  elevation,  7] ,  that  simultaneously  satisfy  the  boundary  conditions. 

Finding  this  simultaneous  solution  is  a  difficult  task.  However,  a  systematic  solution  approach  is  to  iterate 
between  computing  the  flow  beneath  an  approximate  free  surface,  and  using  the  flow  solution  to 
determine  a  better  approximation  to  the  free  surface.  Iteration  continues  until  the  free  surface  for  which 
the  boundary  conditions  are  sufficiently  satisfied  is  determined. 

Now  the  velocity  potential  and  free-surface  elevation  are  decomposed  into  two  parts:  a  current 
approximation  (denoted  with  a  superscript  0),  and  a  refinement  to  this  approximation  (denoted  with  a 

superscript  1).  Specifically,  the  respective  expressions  are0  — >  <f>°  +  <fi'  and  Tf  —¥JJ°  +  if .  Furthermore, 
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the  kinematic  (Equation  2)  and  dynamic  (Equation  3)  boundary  conditions  are  combined  into  a  single 
free-surface  boundary  condition 


E°p\ 


■E°zpl +E1p°z-E'zp°  =E°zp° -E°p 


(6) 


The  use  of  Equation  6  results  in  quadratic  convergence;  however,  up  to  third-order  derivatives  of  the 
velocity-potential  function  are  necessary.  In  addition  to  Equation  6,  a  “radiation  condition”  is  imposed  at 
two  upstream  points  to  keep  waves  from  “radiating”  upstream  and  to  ensure  that  Equation  4  is  satisfied. 
The  particular  radiation  condition  suggested  by  Scullen  [1]  is 


X0xz  +30z  =0. 


(7) 


which  states  that  the  vertical  component  of  the  velocity  decays  in  proportion  to  the  inverse  cube  of 
distance.  By  performing  the  necessary  derivatives  and  expressing  the  results  in  terms  of  <f>°  and  (see 
Appendix  A),  a  system  of  linear  equations  can  be  constructed  whose  solution  determines  the  refinement 
to  the  velocity  potential.  With  this  new  potential,  the  refinement  to  the  free-surface  elevation  is  obtained 
by 


V 


l 


(8) 


From  an  implementation  viewpoint,  the  free  surface  is  discretized  with  a  uniform  distribution  of 
collocation  points,  and  discrete  sources  are  placed  a  fixed  distance  above  each  free-surface  collocation 
point.  Since  Equation  7  is  imposed  in  addition  to  Equation  6  at  two  upstream  collocation  points,  two 
additional  sources  are  used.  One  is  placed  upstream  to  the  first  collocation  point,  and  the  other  is  placed 
downstream  of  the  last  collocation  point  while  maintaining  the  same  x-axis  distribution.  This  can  also  be 
thought  of  as  eliminating  the  first  and  last  collocation  point  that  would  be  directly  beneath  these  two 
sources.  The  number  of  sources  and  their  placement  are  determined  based  on  the  work  and  conclusions  of 
Scullen  [1,2].  While  the  free-surface  collocation  points  are  moved  vertically  according  to  Equation  8,  the 
free-surface  sources  remain  fixed.  Figure  19  shows  a  diagram  of  the  free-surface  representation. 
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Scullen  also  presented  results  for  a  submerged  cylinder  with  prescribed  circulation.  These  simulations 
were  performed  by  specifying  body  collocation  points,  and  placing  sources  85%  of  the  distance  along  the 
radius  to  these  collocation  points.  The  circulation  was  achieved  by  using  a  vortex  placed  at  the  center  of 
the  cylinder  and  its  image  above  the  free  surface.  While  the  reason  for  using  the  vortex  image  is  not 
explicitly  stated,  it  is  believed  that  using  the  image  uniquely  sets  the  upstream  flow  conditions  by  forcing 
the  far  upstream  free-surface  elevation  to  be  zero.  These  simulations  were  reproduced,  and  the 
configuration  and  results  for  one  set  of  the  computations  are  respectively  given  in  Figure  20  and  Figure 


21. 


Figure  20:  Configuration  for  radius-to-depth  ratio  of  0.1333  and  a  Froude  number  of  0.3651. 


Figure  21:  Nonlinear  waves  produced  by  cylinders  with  the  above  configuration  and  varying  circulation. 
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Second,  the  well-established  linearly-varying  vortex-panel  method  [3, 4]  to  compute  the  flow  about  a 
hydro/airfoil  was  implemented  in  a  way  that  is  compatible  with  the  nonlinear  free-surface  solution 
approach  described  above.  The  vortex -panel  method  uses  panels  to  represent  the  surface  of  an  airfoil,  and 
the  distribution  of  vortex  strength  along  each  panel  varies  linearly  from  one  unknown  value  at  one  end  to 
another  unknown  value  at  the  other  end.  As  a  result,  integration  along  the  panel  is  required  to  compute 
each  panel’s  contribution  to  the  global  potential  function  (or  its  derivatives)  at  any  given  point  in  space, 
i.e. 


<t> panel  \r(s)6(x,Z,S)ds  .  (9) 

panel 

Here,  y(s)  represents  the  vortex-strength  distribution  along  the  panel  and  0{x,  z,s)  represents  the  angle 
between  the  line  connecting  the  element  ds  to  the  point  (jc,  z)  and  the  x-axis. 


Each  panel  end  point  is  defined  as  a  point  on  the  surface  of  the  body,  and  the  collocation  points  are  taken 
to  be  the  panel  midpoints.  As  an  implementation  note,  the  N  panels  are  ordered  in  a  clockwise  manner, 
and  the  surface  points  are  obtained  from  an  analytical  correlation  to  the  appropriate  NACA  4-digit  wing- 
section  data.  These  and  the  above-mentioned  concepts  are  illustrated  below  in  Figure  22. 

It  is  easy  to  see  that  the  Neumann  boundary  condition  (Equation  5),  which  is  imposed  at  each  collocation 
point  on  the  airfoil,  involves  first-order  derivatives  of  the  velocity  potential  and  is  linear  with  respect  to 
the  velocity  potential.  However,  the  integral  in  Equation  9  makes  the  derivation  of  the  analytical 
expressions  for  these  quantities  more  complicated.  Typically  the  integration  is  performed  in  a  coordinate 
system  that  has  been  rotated  such  that  the  x-axis  is  aligned  to  the  panel  and  the  first  end  point  of  the  panel, 
(xpi,  zpi),  is  at  the  origin.  This  rotated  coordinate  system  will  be  referred  to  as  the  panel  coordinate  system 
(PCS).  Given  the  angle,  J3 ,  between  the  panel  and  the  x-axis  in  the  global  coordinate  system  (GCS),  this 
is  achieved  through  applying  the  following  transformation  to  the  x-  and  z-coordinates  of  the  point  of 
interest: 


xpcs  =  (xgcs  ~  xPi  )Cos(P)  +  (zGCS  -  zpl  )Sin{P)  and  (10) 

Zpcs  =  ~{XGCS  ~  Xp\  +  ( Z-GCS  —  Zp\  )cVw(/?)  (11) 
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The  potential  can  then  be  obtained  by  simply  performing  the  integration  and  arranging  the  result  in  terms 
of  the  unknown  end-point  vortex  strengths.  The  derivatives  of  the  velocity  potential,  with  respect  to  the 
GCS,  can  also  be  obtained  through  using  Equations  10  and  1 1  and  applying  the  chain  rule  of 
differentiation.  (Note:  since  s  is  a  “dummy”  variable,  differentiation  of  Equation  9  with  respect  to  the  PCS 
can  be  simply  taken  inside  of  the  integral  and  is  only  applicable  to  0(x,  z,  s ) .)  Appendix  A  provides 

explicit  expressions  for  the  first-order  derivatives  of  the  vortex-panel  velocity  potential  in  terms  of  the 
end-point  vortex  strengths. 

Since  the  boundary  condition  on  the  surface  of  a  hydrofoil  is  linear  with  respect  to  the  velocity  potential, 
the  flow  solution  is  typically  obtained  by  solving  a  single  linear  system  of  equations  rather  than  iterating. 
To  reconcile  this  difference  in  approaches,  the  velocity  potential  corresponding  to  a  vortex  panel  was  also 
decomposed  into  an  approximation  and  a  refinement  to  the  approximation.  With  this  in  place  the 
remaining  task  was  to  implement  a  Kutta  condition  to  close  the  system  of  equations  and  to  provide  a 
means  to  establish  the  amount  of  circulation  required  to  properly  compute  the  flow  about  a  lifting  body. 
Typical  Kutta  condition  implementations  force  the  sum  of  the  trailing-edge  vortex  strengths  to  be  zero, 
but  this  will  not  take  into  consideration  any  other  sources  of  velocity  potential  that  might  be  present.  So, 
an  alternative  approach  [6]  was  chosen  which  requires  the  flow  at  a  point  just  behind  the  trailing  edge  of 
the  hydrofoil  to  be  tangent  to  the  extended  symmetry  line  of  the  airfoil. 


Figure  23:  Comparison  of  flow  solutions  around  a  NACA  0012  with  an  altered  and  unaltered  panel 

method  using  200  panels. 


To  verify  the  modified  method’s  implementation,  solutions  (see  Figure  23)  for  the  flow  around  a  NACA 
0012  airfoil  at  an  angle  of  attack  of  5°  were  computed  using  both  an  unaltered  vortex-panel  method  code 
from  FDLIB  [7]  and  the  described  modified  method.  Both  simulations  used  200  panels  to  represent  the 
airfoil.  The  pressure  coefficient  distributions  shown  in  Figure  23  are  in  very  good  agreement.  Both 
methods  display  undesirable  characteristics  near  the  trailing  edge.  Wauquiez  also  points  this  out  and 
describes  another  type  of  Kutta  condition  using  a  wake  panel  [4]  that  improves  the  solution  near  the 
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trailing  edge.  This  Kutta  condition  implementation,  however,  was  not  incorporated  into  this  work.  Also 
note  that  the  analytic  correlation  that  was  used  for  the  airfoil  does  not  close  at  the  trailing  edge;  this  could 
also  explain  some  undesirable  behavior  near  the  trailing  edge. 

Third,  the  two  methods  were  combined  through  completing  several  tasks.  Expressing  Equation  6  in  terms 
of  the  velocity  potential  (See  Appendix  A.)  reveals  that  up  through  third-order  derivatives  are  required  for 
the  nonlinear  free-surface  method.  Only  first-order  derivatives  are  necessary  for  the  panel  method. 
Therefore,  the  first  task  was  to  derive  the  necessary  higher-order  derivatives  of  the  velocity  potential  for  a 
panel.  This  was  accomplished  by  simply  following  the  same  approach  that  was  used  to  derive  the  first- 
order  derivatives,  and  Appendix  A  provides  explicit  expressions  for  the  necessary  derivatives.  Another 
task  necessary  to  combine  the  methods  was  to  mirror  the  vortex  panels  above  the  free  surface.  This  is 
directly  analogous  to  mirroring  the  single  vortex  used  in  the  simulations  of  a  cylinder  with  circulation, 
except  now  we  are  dealing  with  vortex  panels.  The  last  task  was  to  incorporate  a  method  to  simulate  a 
finite-depth  configuration.  This  could  be  achieved  by  placing  sources  of  unknown  strength  outside  of  the 
domain  and  below  the  fluid  bottom.  The  unknown  source  strengths  could  then  be  computed  to  satisfy  the 
Neumann  boundary  condition  along  the  fluid  bottom.  A  better  approach  that  does  not  increase  the  size  of 
the  linear  system  of  equations  to  solve  is  to  simply  mirror  every  singularity  about  the  fluid  bottom.  The 
latter  approach  implicitly  satisfies  the  Neumann  boundary  condition  along  the  fluid  bottom,  and  it  was 
utilized  in  this  work. 

Results 

The  method  described  above  was  utilized  to  perform  simulations  corresponding  to  configurations  used  in 
three  sets  of  experiments.  The  first  set  was  that  of  Duncan  [5];  the  second  was  performed  by  Furey,  et  al 
[8];  and  the  third  by  Walker,  et  al  [9].  Due  to  limitations  of  the  method  to  produce  solutions  for  large 
wave  amplitudes,  results  were  obtained  for  only  some  of  the  experiments  within  the  first  two  sets  of 
experiments.  The  simulations  that  produced  converged  solutions  typically  involved  only  non-breaking 
waves  or  at  most  mildly  breaking.  The  simulations  that  produced  converged  results  are  described  and 
presented  below. 

Duncan  [5]  performed  experiments  where  a  hydrofoil  with  a  NACA  0012  wing  section  was  moved 
through  the  water.  The  hydrofoil  with  a  chord  of  0.203  m  was  given  a  5°  angle  of  attack  with  the  pivot 
point  at  mid  chord.  The  speed  of  the  hydrofoil  was  0.8  m/s,  and  the  hydrofoil  was  placed  0.175  m  above 
the  fluid  bottom.  In  addition  to  these  parameters,  experiments  were  performed  with  differing  free-surface 
heights  above  the  hydrofoil.  The  free-surface  heights  studied  include:  0.261  m;  0.236  m;  0.210  m;  0.193 
m;  and  0.185  m.  As  the  hydrofoil  approaches  the  free  surface  the  resulting  wave  amplitude  increases  and 
wave  breaking  occurs.  Converged  solutions  were  obtained  for  the  first  three  free-surface  heights.  For  each 
solution  two  plots  are  presented  below.  The  first  plot  shows  the  free-surface  elevation  (with  the 
experimental  data,  if  available),  and  the  second  shows  the  pressure  coefficient  distribution  on  the 
hydrofoil  compared  with  that  of  a  hydrofoil  in  an  infinite  domain.  Figure  24  through  Figure  29 
respectively  present  the  results  for  the  first  three  free-surface  elevations  above  the  hydrofoil. 
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Figure  24:  Free-surface  elevation  for  Duncan’s  depth  =  0.261  m  case. 
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Figure  25:  Negative  pressure  coefficient  distribution  for  the  above  case. 
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Figure  26:  Free-surface  elevation  for  Duncan’s  depth  =  0.236  m  case. 
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Figure  27:  Negative  pressure  coefficient  distribution  for  the  above  case. 
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Figure  28:  Free-surface  elevation  for  Duncan’s  depth  =  0.210  m  case. 


Figure  29:  Negative  pressure  coefficient  distribution  for  the  above  case. 


The  four  experiments  presented  by  Furey  et  al  [8]  also  utilize  a  hydrofoil  with  a  NACA  0012  wing 
section,  but  its  pivot  point  is  at  the  quarter-chord  position  and  the  flow  moves  past  the  hydrofoil.  Here  the 
fluid  bottom  is  fixed  at  0.9  m  below  the  hydrofoil  which  has  a  chord  of  1 .8  m  and  an  angle  of  attack  equal 
to  5.5°.  Both  the  undisturbed  free-surface  elevation  and  the  free-stream  velocity  are  allowed  to  vary 
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between  experiments.  There  seems  to  be  some  contention  regarding  whether  the  reported  flow  velocities 
are  correct;  velocities  obtained  from  an  acoustic  Doppler  velocimiter  (ADV)  were  also  available.  The 
ADV  velocities  were  considerably  larger.  While  converged  solutions  were  obtained  using  both  velocities 
for  the  first  experiment,  a  solution  could  be  produced  for  only  the  lower,  non-ADV,  velocity  for  the 
second  experiment.  Figure  30  and  Figure  3 1  present  results  for  the  first  experiment  using  each  free-stream 
velocity,  1.86  m/s  and  2.10  m/s.  and  an  undisturbed  free-surface  elevation  of  1.79  m  above  the  hydrofoil. 


Figure  30:  Free-surface  elevations  for  the  first  case  of  Furey  et  al  using  two  different  velocities. 


Figure  31:  Negative  pressure  coefficient  distribution  for  the  above  case. 
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Nonlinear  Free-Surface  Expressions 


This  section  presents  the  components  of  the  free-surface  boundary  condition  in  terms  of  the  potentials  <p° 
and  (f)x .  It  is  assumed  that  the  pressure  above  atmospheric  is  defined  as 

+WJ  -ul)+gz 

P  2 

where  p  is  the  fluid  density,  g  is  gravity,  and  U„  is  the  free-stream  velocity.  It  is  also  assumed  that  E  is 
defined  as  the  total  derivative  of  the  pressure,  i.e. 

e = {<p*  y  + 2 ^z^xz + (0z  y  ^zz + • 

With  these  definitions  in  place  and  remembering  that  z  =  TJ°  on  the  free  surface,  the  expressions  for  the 
free-surface  boundary  condition  can  be  derived.  While  they  are  presented  by  Scullen  [2],  they  are 
presented  here  as  well  for  completeness.  The  various  components  with  the  corresponding  expressions 
used  in  this  work  follow. 

P°  -U1-  )+sn° 

p"  =  «*,Vi+«>X+* 

E'  =  +  fJt  V;  +  fe?  1  €  +  WWn  + 

2 +|*)rt!+ 

k = (/Jc + wyji + wwa + wyji + 


K = 2  tec + M  k + (€  1 C + 
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2k?+*i& +tef+w«k +(g+wxk+tef€. 


First-Order  Vortex-Panel  Velocity-Potential  Derivatives 

As  was  stated  above,  the  first-order  derivatives  of  the  velocity  potential  for  a  vortex  panel  with  linearly- 
varying  vortex  strength  are  obtained  by  transforming  to  a  coordinate  system  relative  to  the  panel,  denoted 
PCS,  and  utilizing  the  chain  rule  for  differentiation  to  compute  the  derivatives  in  the  global  coordinate 
system  (GCS).  This  section  provides  the  results  of  such  a  procedure. 

Recalling  the  definition  of  the  transformation  from  the  GCS  to  the  PCS  described  by  Equations  10  and 
1 1 ,  one  can  obtain  the  following  expressions  for  the  derivative  of  the  PCS  with  respect  to  the  GCS  as 


5* 

dx 


PCS 


GCS 


=  Cos(/J) ,  ^ 


dz, 


GCS 


=  Sin(j5 ) 


dZpcs 

U^GCS 


=  -Sin(j3) ,  and 


^  =  Cos{p). 

vZgcs 


Now,  the  derivative  of  the  velocity  potential  with  respect  to  x  in  the  GCS  becomes 


dx 


d0_  =  d(p  dXpcs  +  df  dzpcs  =  Cos(^Jl - Sin{py 


GCS 


dxpcs  dxGCS 


dz 


pcs  uagcs 


dx 


PCS 


dz 


PCS 


and  the  derivative  with  respect  to  z  in  the  GCS  is 


d(j>  _  d(f>  dxpcs  d(f>  dzPCS 
dz gcs  dxPCS  dzgcs  dzPCS  dzGCS 


=  Sin(P)-^-  +  Cos(P)y4-. 

dXpcs  dz,  pcs 


The  remaining  task  is  to  express  the  derivatives  of  the  panel’s  velocity  potential  in  terms  of  the  end-point 
vortex  strengths,  y.  and  y  .+1 .  In  the  PCS,  the  expression  for  6(x,z,s)  becomes 


f 

6{x,  z,  s)  =  ArcTan 


\ 


46 


DW200392.8.F.Nov04 


W  _  r, 


3 ^  pcs  2  j dc  2 


f 

f  \ 

\ 

Zpcs 

Zk 

+  ( XPCS  ”  X2  )(^2  “  ^1  ) 

V 

KV\  ) 

) 

?0+. 


2mc, 


3f»  _  J'/ 


3^  pcs  2m2 


r 

(  r  \ 

Zpcs 

r2 

+  x  res  {@2  ~  0\ ) 

V 

1  J 

\ 

XPCS  ~  x 

2)ln 

—  \-ZpCs{0 2~ 

U  J 

r** 


2nx , 


f 

r  \ 

\ 

XPCS 

Zk 

“  Zpcs  {@2  “  @\  )  +  X2 

V 

Kr\  J 
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The  following  diagram  defines  the  newly  introduced  terms. 

Higher-Order  Vortex-Panel  Velocity-Potential  Derivatives 

In  a  similar  fashion  the  higher-order  derivatives  of  the  vortex-panel  velocity  potential  can  be  derived. 
Once  again,  the  coordinate  transformation  definition  is  used  along  with  repeated  use  of  the  chain  rule. 
Further  simplifications  are  also  incorporated  by  using  the  facts  that  the  velocity  potential  satisfies  the 
Laplace  equation,  i.e.  V2<p  =  0,  and  that  the  orders  of  differentiation  are  interchangeable. 

With  all  of  this  in  mind,  expressions  for  the  necessary  higher-order  derivatives  are  given  below: 

=  [Cos\P)  -  Sin\p ihzr-%: - 2  Sin(P)Cos(p)  — ^ - , 

ox  pcs  OZ  pcs 


^XGCS^XGCS 


7)jr 

UAPCSUAPCS 


- =  2 Sin(P)Cos(p)~  +  (cos2  (P)  -  Sin2  (P))—^— , 


GCS^GCS 


dxpcs^Xpcs 


dxPCndz 


pcs  PCS 


- =  (Cos3  (P)  -  3  Cos(P)Sin2  (/?))- - - + 


^XGCS  ^ XGCS  ^ XGCS 


d xpcs  ^x  pcs  dxpcs 


(: iSin(p)Cos2(P )  -  Sin3(p))~ - - ,  and 

OZ  pcs  OZ  pcs  OZ  pcs 

- - - =  (sin3  ip)  -  3  Sin(P)Cos2  (/?))- - -f*- - + 

OZ  CCS  UZGCS  ®ZGCS  °XPCS  °XPCS  °xpcs 

(Cos3{p)  -  3 Cos(p)Sin2(P))~ - - . 

OZ  pcs  OZ  pcs  OZ  pcs 


As  before,  the  final  task  is  to  compute  the  panel  velocity-potential  derivatives.  With  the  same  definitions 
as  above,  their  expressions  follow: 
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